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3-D MU LTILATERATION : A PRECISION GEODETIC MEASUREMENT SYSTEM 


ABSTRACT 

The systems analysis and laboratory demonstration described in this 
report indicate that a new technique of satellite geodesy, to be called 3-D Mul- 
tilateration, can determine the relative three dimensional coordinates of 
ground stations within 1 centimeter over baselines of 20 to 10, 000 kilometers. 
With this high accuracy, several crucial geodetic applications become possible. 
These applications include: earthquake hazards assessment, precision survey- 

ing, plate tectonics, and orbital applications. 

Achievement of such accuracy can be attained through use of pulsed lasers 
to obtain simultaneous slant-ranges between an ensemble of ground stations and 
a moving r etror eflector whose trajectory is known a-priori only to the accuracy 
necessary for aiming the lasers. Specifically, the positions of the satellite or 
airplane carried retr oreflector s are eliminated from the equations which govern 
determination of station locations. However, once the station locations are 
determined, the trajectory of the r etroreflector can be obtained as a direct 
by-product. 

Numerical analysis has shown that suitably chosen multistation configura- 
tions result in well-conditioned solutions, with very small error magnification 
of the inherent ranging errors occasioned by the hardware subsystem. 

Laboratory tests have demonstrated that a laser hardware subsystem with 
a ranging accuracy of 3 centimeters can be built from commercially available 
components. By 1975, at the latest, an accuracy level of 1 centimeter can be 
achieved. Costs of the proposed systems are lower than other proposed sys- 
tems on a station construction and implementation basis. 
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SECTION I 


SUMMARY 


2. This Section briefly summarizes the basic concepts and investigations 
of the 3-D Multilate ration system study performed at JPL under contract 
number NASA 7-100. The Summary Section is intended to be a detachable unit 
from the main body of this report. 

1. 1 APPLICATIONS OF THE 3-D MULTILATERATION TECHNIQUE 

Earthquakes are the result of sudden slips between two blocks of the 
Earth's crust. The entire Earth is covered by large blocks which may be in 
steady motion relative to each other at points distant from their boundaries. 

At the boundaries, i. e. , at seismic zones, the blocks are held together by 
friction. When the driving stresses exceed the breaking strength of the fric- 
tional contact, abrupt motion ensues and an earthquake occurs. However, the 
Earth is a deformable body and in the earthquake case, it deforms non- 
elastically in creep. Creep strains build up in regions astride earthquake 
faults, and are relieved, at least in part, by the occurrence of an earthquake 
on the fault. Occasionally, some creep is observed on the faults themselves. 

A knowledge of the strain configuration in the Earth's crust, both at a fault and 
at points distant from the fault, is an essential ingredient to a program of 
learning how to estimate the location and the epoch of major earthquakes. 

There are no specific problems in the measurement of creep along faults, 
since the motion takes place over relatively short distances. To evaluate the 
strain accumulation in a large region surrounding a fault, it is necessary to 
measure precisely the relative motion of points in various parts of this region. 
If this motion is to be evaluated within a period of one to two years, it will be 
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necessary to measure distances as large as 10, 000 km with an accuracy of 3 cm 
or better. Until this moment there has been no technique capable of making 
such measurements. 

In this report a new technique, 3-D Multilateration, will be described. 

This technique appears to provide not only a one centimeter distance measure- 
ment accuracy capability, but also determines each station position in three 
dimensions and thereby allows evaluation of the different components of strain 
buildup. Furthermore, since the 3-D Multilateration technique will require only 
hours'to obtain raw data and process this data into geophysically significant 
information, it may be possible to detect the presence of rapidly occuring 
ground creep which often gives a premonitory indication of earthquakes. 

It should be noted that the 3-D Multilateration technique can be used in 
geophysical analysis not directly related to earthquake prediction. This tech- 
nique is sufficiently flexible to permit a wide variety of applications, as illus- 
trated in Figure 1-1. 



Figure 1-1. Applications of 3-D Multilateration 
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In summary, if the performance of the 3-D Multilateration technique 
demonstrated in this report is realized, not only will science be provided with 
a valuable tool for investigation of tectonic phenomena, but there is a significant 
pos sibility that geophysicists will be able to predict the location, and perhaps 
the approximate time, of future earthquakes. 

1. 2 DESCRIPTION OF THE 3-D MULTILATERATION TECHNIQUE 

The 3-D Multilateration technique can be implemented through use of a 
number of ground stations which nearly simultaneously transmit laser pulses to 
reflectors on a moving vehicle, e. g. , an airplane or a satellite. Each station 
evaluates station-to- vehicle range by measuring the time interval between 
transmission of the original pulse and reception of the reflected pulse. Simul- 
taneous range measurements are then processed to yield relative station 
locations in 3 dimensions. 

The components which are necessary to define the Precision Geodetic 
Measurement System are illustrated in Figure 1-2. 



Figure 1-2. Two Basic Components of Overall System 


See Para. 1. 5 for definition of the phrase "relative station locations. " 
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An important advantage of multilateration techniques is that no information 
concerning the vehicle trajectory is required to obtain station coordinates rela- 
tive to one another. Thus, uncertainty of the satellite ephemeris has no effect 
on the deduced station coordinates. Indeed, it is possible to determine the tra- 
jectory of the vehicle in the coordinate frame of the stations to the centimeter 
level. In summary two kinds of information are provided by 3-D Multilateration: 

• The three-dimensional coordinates of all stations utilized in the 
operational configuration, and 

• The precise positions of the moving r etroreflector mounted on the 
vehicle. 

1. 3 THREE SYSTEMS USING MULTILATERATION TECHNIQUES 

Three separate systems, which are possible, will be discussed in this 
report. These systems will require the use of six, four and three stations for 
operational implementation. Specifically it is demonstrated that solutions 
for the station coordinates are possible if: 

• Six stations simultaneously make 6 range measurements (between 
station and vehicle) at four independent times, 

• Four stations simultaneously make 4 range measurements at six 
independent times, 

• Three stations (which are aligned) simultaneously make 3 range 
measurements at two independent times. 

The above three solutions give rise to what will be referred to in this report as 

the 6 station system, the 4 station system, and the 3 station system. 

The process of making simultaneous range measurements from any num- 
ber of stations to a vehicle, i. e. , to an airplane or satellite (at a single instant 
of time) is called a strike . Hence, to obtain a solution for the station coordin- 
ates what is required is a sequence (timewise) of strikes. The rules described 
above are stated as the minimum requirement for a deterministic solution to 
the 3-D Multilateration problem. A much longer sequence of strikes would be 
used for purposes of statistical improvement Hence in more detail there exists: 
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The Six Station Solution (Figure 1-3) which states that: 


Six ground stations which make a sequence of four or more simultaneous 
strikes (which are separated in time by any arbitrary duration) to one or two vehi- 
cles, will give rise to a mathematically and numerically stable solution for the 
determination of relative three-dimensional station locations. The station coor- 
dinates can be obtained without using the vehicle coordinates within the defining 
equations, and without any restriction on the distance between the stations. The 
vehicle position and velocity can be obtained as a by-product of the calculation. 



See Para. 1. 5 for the definition of relative station locations. 
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The Four Station Solution (Figure 1-4) which states that: 


A solution for the relative station locations can be obtained by using four 
ground stations which make a sequence of six or more simultaneous strikes 
(which are separated in time by any arbitrary duration) to one or more vehicles, 
if the fourth station is located out of the plane of the other three stations by a 
reasonably large distance, such as might be obtained in establishing interconti- 
nental reference locations. As with the previous solution the position of the vehi- 
cle does not enter into the calculations, but is in fact obtained as a by-product. 



Figure 1-4. Four Station Solution 
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• The Three Station Solution (Figure 1-5) which states that: 

The relative distance between two ground stations can be obtained by 
aligning a third station in a straight line with respect to the first two stations 
and taking a sequence of two or more simultaneous strikes (which are separated 
in time over any arbitrary duration) to a vehicle. Again, no satellite ephemeris 
information is needed in the calculation of station locations. 



ONLY 

BASELINE 

DISTANCES 

BETWEEN 

STATIONS 

DETERMINED 


Figure 1-5. Three Station Solution 


JPL Technical Memorandum 33-605 


1-7 





In closing this Section it should be mentioned that other system 
combinations are possible. For example, in the four station solution it is 
possible to fix two of the baselines between the stations and arrive at a suitable 
solution for nearly coplanar configurations. In fact, the four station system 
previously described, can then be used for local and intercontinental applica- 
tions. The baselines can be obtained by several methods which will be described 
shortly. Hence, though not much emphasis is placed upon these systems there 
also exists a myriad of: 

• Hybrid Multilateration Solutions. 

These and other solutions can be obtained using range data (as previously 
described) or by use of range difference data, i. e. , the difference in-range 
between the vehicle and two stations. 

1.4 HARDWARE SYSTEM AND ACCURACY DEMONSTRATION 

The hardware subsystem for measuring station-to-vehicle ranges is iden- 
tical for all station configurations. In order to obtain high accuracy range 
measurements, the subsystem utilizes a new type of pulsed laser. This laser is 
a mode-locked, Q- spoiled ruby laser, and has the capability of emitting very 
short duration (0. 1 ns or less) light pulses. The subsystem employs a tracking 
mount to aim transmitting and receiving telescopes at the vehicle. A measure- 
ment is made of the duration required for a laser pulse to make a round trip 
from the station to the vehicle borne retroreflector . This duration, measured 
with a resolution of 0. 1 ns, is used in conjunction with an atmospheric model to 
calculate the range to the vehicle. 

Each ground station contains an X-Y tracking mount for steering the two 
telescopes, a small computer to direct the tracking mount towards the satellite, 
timing circuitry, recording equipment and power supplies. 
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In order to satisfy the requirement for "simultaneous" ranging, the clocks 
at each station must be synchronized to 3 psec. Such synchronization is easily 
achievable using low-cost components. Synchronization of laser firing to 1 ms 
accuracy is adequate and well within the state of the art; variations of laser 
firing within this range are compensated for by time-tagging range measure- 
ments, and interpolating between successive measurements so as to obtain 
effective simultaneity among stations. 

In order to evaluate the measurement errors produced by current hard- 
ware, a demonstration ranging subsystem was assembled using commercially 
available components. The hardware configuration simulated a long distance 
ranging system over short path lengths by using attenuated return signals. The 
results obtained with this system show that hardware related errors in range 
measurement can be made acceptably small, typically less than 2 cm. Future 
systems, fabricated circa 1975, can be expected to operate with 1 cm ranging 
accuracy for satellites placed in orbits yielding nominal system efficiency. 

1. 5 ANALYSIS OF SYSTEM ERRORS 

It is emphasized that the system errors caused by the satellite. Earth con- 
stants, and orbital perturbations do not enter into the multilate ration process. 

In fact, since the proposed techniques are independent of the location of the 
retroreflector, the only error sources which enter into these techniques are as 
follows : 

• Multiple reflection error, 

• Bias error due to atmospheric delay, 

• Random error due to atmospheric turbulence, 

• Random equipment measurement error, 

• Equipment bias error. 

It is shown in this study that a ranging subsystem can be fabricated in which the 
net effect of these errors will lead to a ranging accuracy of 1 cm. 
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The geometric coordinate system used in these studies is a relative coor- 
dinate system, in which Station 1 is placed at the origin, arbitrary placement of 
Station 2 fixes the X axis, and arbitrary location of Station 3 defines the X-Y 
plane of the coordinate system (see Figure 1-6). 

This coordinate system can be linked to the inertial geocentric coordinate 
system only via the geographic coordinates of the first adopted station. The 
standard deviations in the latitude, longitude and elevation of Station 1, estab- 
lished by other methods distinct from multilateration, are required inputs 
to the process of computing errors in the inertial coordinate system. 

These error sources, however, are not an integral part of the proposed 
geometric laser ranging technique which has been studied at the Jet Propulsion 
Laboratory for the determination of relative station locations. 

If the method of multilateration is used to determine the trajectory within 

the standard inertial system, then the following errors must be considered: 

• The geographic coordinates of the first station, 

• The astrodynamic constants, 

• Timing errors. 


z 
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Specifically, these errors will be transferred only to the coordinates of the orbit 
in the inertial frame. Hence only the coordinates of the orbit will be corrupted 
by the analytic process but not the relative station locations.* If the trajectory 
of the orbit is determined in the geometric frame these errors will, of course, 
vanish. 

1. 6 ANTICIPATED PERFORMANCE OF THE PROPOSED SYSTEMS 

The comparison of the three proposed systems requires examination of 

three basic factors: 

• System Cost, 

• Operational Cost, 

• Accuracy. 

1. 6. 1 Cost Estimates 

Station fabrication costs will be moderate, with each station expected to 
cost approximately $275, 000 when fabricated in quantities of six. 

While operational costs have not been estimated in this report, automatic 
operation of the hardware subsystem has been included in the preliminary design. 

The software subsystem costs are a small fraction of the overall cost 
due to the inherent simplicity of 3-D Multilate ration technique. An operational 
software system to process the raw data can be developed for $150, 000. 

1.6.2 Accuracy 

It is shown in this report that the three and four station systems have 
certain limitations. Specifically the three station system will function only 
if the three stations are collinear. Furthermore, this system determines inter- 
station baselines only. On the other hand, multilate ration systems using four 
or five stations do not function for local separation between the stations. 
Specifically, the equations defining the three dimensional coordinates become 

*The laboratory units of distance discussed in this report are ultimately 
dependent on the adopted value for the velocity of light. 
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singular unless the stations are placed at intercontinental distances. However, 
it is proved that if six or more stations are used to define the system, the inher- 
ent mathematical singularity vanishes. This important fact has the effect of 
increasing the accuracy of the station location determinations whether the six 
station system is used in a local or global capacity. For observation sets of 
100, 1-cm accuracy measurements, the performance of each proposed system 
is summarized below. 

• For the six station system, the errors occasioned upon the three 
coordinates of each station are of the order of 1 centimeter for well- 
scattered local and/or global station locations. 

• For the four station system, using a wide baseline separation 
(3000-5000 km) between the stations, an accuracy of the order of 
5 to 10 centimeters can be obtained. 

• If three stations in a straight line are used to determine the 
distance between the stations, e.g. , in the three- station system, 
an accuracy of 3 centimeters can be obtained in determining the 
baseline between the first two stations. 

1.7 OPERATIONAL CONSIDERATIONS 

As will be explained, the proposed systems require the use of an elliptic 
orbit, or two circular orbits. These requirements are not a barrier and are 
not expensive. The elliptic orbit approach is recommended as the simplest, 
e.g. , a 400 x 1000 km orbit with inclinations between 45° and 60° is acceptable 
for stations with 100 to 1000 km separation. 

The proposed system outlined herein can be implemented by using a 
octahedral array of reflectors to achieve a strict point source reflection. The 
array can be attached to any available satellite whose overall mission objective 
might be considerably different from the experiments proposed herein. 
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Atmospheric attenuation is not expected to significantly degrade system 
performance. Other work (See Section 8) indicates that atmospheric effects 
can be corrected to a residual error of 2 mm at the zenith, and approximately 
1 cm at 70° from the zenith. 

Since each station uses small telescopes and a simple tracking mount, 
the stations can be made transportable and can be moved to various locations. 
Thus, a single set of stations can be used to map a large number of locations. 

It has been already mentioned that synchronization of the laser firing to 
1 ms is adequate, and that data-interpolation can be obtained from polynomials. 
The easily met constraint of the station clock synchronization to 3 psec, which 
is necessary for time-tagging the observations, presents no problem. 

In real time operations the station up-times due to common visibility 
durations as a function of weather must be considered. Safety considerations, 
as with all laser systems, must be taken into account; e. g. , it has been shown 
that laser damage to the unshielded human retina can result if a direct laser 
strike is sustained by a pilot whose retina crosses the laser beam. 

1.8 COMPARISON WITH OTHER TECHNIQUES 

There are three other techniques which are presently in existence for 
precision geodetic surveying. They are: 

• The Geodolite 

• The Goddard Laser Ranging System 

• Very Long Baseline Interferometry 

Spectra-Physics of Mountain View, California manufactures a laser- 
ranging device known as the Geodolite which is capable of measuring the distance 
between points on the Earth's surface. The Geodolite can measure only 


JPL Technical Memorandum 33-605 


1-15 



line-of-sight distances, but provides excellent accuracy over its severely 
constrained range. For example, the U.S. Geological Survey has used a 
Geodolite to measure baselines of 30 km with accuracies of 1 cm or better; 
however, in order to achieve this accuracy, the USGS combined a number of 
partially redundant measurements from different benchmarks, and then cor- 
rected these raw measurements with temperature and humidity data obtained 
by overflying the baselines with a helicopter. In summary, although the 
Geodolite appears to be an excellent device for obtaining high accuracy distance 
measurements over short baselines, it is not an acceptable tool for measuring 
the 3 -dimensional components of distance over long baselines. Consequently, 
the Geodolite is suitable for investigations similar to those anticipated for our 
3 -station system, but is completely inapplicable to the long distance, 

3 -dimensional measurement problems capable of being solved by 6- station 
s ystem. 

The second comparable technique is that developed by Goddard Space 
Flight Center. The GSFC system also employs laser ranging to a retroreflecting 
satellite, but requires only two ground stations. Three-dimensional station 
locations are obtained by combining the station-to- satellite range measure- 
ments with a highly precise model of the force field acting upon the satellite. 
Goddard presently obtains ranging measurements with a standard deviation of 
approximately 35 cm. The resultant residuals in 3 -dimensional station locations 
have standard deviations of about 50 cm, and the overall accuracy of the station 
location is about 2 meters. GSFC anticipates improving the accuracy of its 
laser ranging equipment so as to yield ranging errors of 1-3 cm (similar to the 

❖Although it would be conceptually possible to augment the Geodolite measure- 
ments with measurements of vertical motion provided by tiltmeters or gra- 
vimeters, such systems are not feasible in practice (See Section 3). 
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ranging errors predicted for the proposed systems). However, the Goddard 
system is extremely sensitive to inaccuracies in modelling the force field 
acting on the satellite, and it is doubtful if such modelling can be improved 
sufficiently to permit evaluation of 3 -dimensional station locations with accur- 
acies of 10 cm or better in the foreseeable future. 

Very Long Baseline Interferometry is a method which by means of 
measuring radio signals from extragalactic radio sources can deduce station 
locations. This method has promise. Accuracies of the order of 10 cm are 
possible. In the 1975 time frame the limit of accuracy will probably approach 
5 cm. The disadvantage of this method is its dependence on state of the art 
frequency standards and digital recording systems and a reliance on major radio 
tracking facilities. This implies increased cost in real time operations. Fur- 
thermore, unlike the 3-D laser ranging system proposed in this report, the 
VLBI system relies heavily on data processing. The software expenses are at 
least an order of magnitude larger than those incurred in the 3-D Multilatera- 
tion process. 

However, when VLBI is used to obtain range differences (see Para- 
graph 1. 3) from an orbiting satellite, then the use of multilate ration techniques 
will drastically alter the accuracy and software aspects of VLBI in a positive 
sense . * 


*The application of range difference multilate ration is currently under 
preliminary investigation. 
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1.9 


CONCLUSIONS 


A precision geodetic measurement system which enables three 
dimensional station coordinates to be determined relative to each other within 
a few centimeters has been demonstrated as being feasible. 

It is concluded that a system using lasers and six operational stations 
has the potential of resolving the coordinates of stations to within 1 cm. 

Having established a system which has such powerful resolution poten- 
tial, serious work can now begin on the development of a workable stress/strain 
geodetic model. This model would have profound influence in the area of earth- 
quake hazards assessment. 

The cost of each laser subsystem has been estimated to be $275, 000. 
The software development costs for the operational software package that will 
process the data are inexpensive. 

The results of the system analysis performed herein indicate that no 
pertinent barriers against implementation of such a system exist. In fact since 
orbits have already been trilaterated using ordinary tracking equipment [l.l], 
the six station concept could be tested at this date by using an additional three 
radio tracking stations. 

Furthermore many other proposed systems, each with unique advantages, 
i. e. , hybrid methods, are also available to aid in geodetic applications. 
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SECTION II 


HISTORY OF PROJECT 

The use of geometric principles to determine the ephemeris of a satellite 
and other related parameters is not a new idea. However, it appears that most 
geometric techniques were abandoned in favor of dynamical methods due to the 
natural evolution of the artificial earth satellite. 

2. 1 BACKGROUND 

Late in 1957, the first man-made satellite was lofted into a low Earth orbit. 
The limited view periods of low satellites as observed by Earth ground stations 
imposed an important operational constraint, namely: the smaller the number 
of ground stations, the lower the overall system cost. Operating under this 
constraint, the idea of using simultaneous geometric ranging from at least 
three stations in order to trilaterate, that is, to geometrically determine three 
triangles whose subsequent solution would yield the position of the satellite, was 
abandoned. Indeed, the dynamic methods of orbit determination were more 
suited to handling the pressing problem of satellite ephemeris propagation at 
those early times. 

As time passed, satellites were lofted into higher orbits and finally into 
the practical and useful geosynchronous orbits now so common for satellite com- 
munication systems. The geosynchronous equatorial satellite, as is well known, 
would appear to hang motionless over a given equatorial station. Since ground 
stations placed on baselines of an intercontinental dimension would have common 
visibility with respect to each other, this advantage was quickly exploited to 
yield a communications relay system. Furthermore, at that date, since station 
nets were well dispersed throughout the globe, visibility of the geosynchronous 
satellites by at least three stations became a common occurrence. This 
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advantage was recognized and the geometric method of trilateration came into 
favor again. Operational software to perform the orbit determination was 
developed and tested by assuming that the location of the stations relative to the 
geodetic datum was known. Analysis indicated that the errors in the satellite 
positions and velocities obtained via trilateration would be smaller than the 
estimates provided via dynamical theory. This was due to the fact that all the 
dynamical constants and their corresponding uncertainties were not required in 
the geometric modeling process. 

2.2 RECENT DEVELOPMENTS 

One of the pressing problems which still remained unanswered until 
recently was concerned with the manner in which station locations themselves 
would be determined. In essence, the question arose of how both the 

position and velocity of the satellite and the three-dimensional station 
locations could be computed simultaneously. 

Early in 1971 some preliminary investigation of this problem was per- 
formed at the Jet Propulsion Laboratory because it appeared that by using laser 
ranging equipment the data, i. e. , the ranges between the station and the satel- 
lite, could be determined to a high degree of accuracy. The objective was to 
estimate Earth station positions, at least the relative positions, an order of 
magnitude more accurately than was presently possible for purposes of geodetic 
applications. 

The results of that investigation showed that if simultaneous ranging is 
performed from a group of four stations at six independent times, then, in a 
specially chosen coordinate system, the positions of the four stations and the 
six trajectory points could be theoretically determined. A second indenture 
analysis discovered a serious problem: namely the station ground separation 
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needed to be on an intercontinental scale, because if all four stations were 
contained in a plane, the mathematical solution becomes singular, i. e. , 
a solution can not be found. Specifically, the geodetically important situation 
of stations with short baselines in a nearly coplanar configuration, e. g. , such 
as in the region surrounding the San Andreas fault in California, could not be 
implemented by the four station method. 

Subsequent investigations discovered the interesting fact that if two mea- 
surements of the range are taken from each of three stations which are all in a 
straight line, then the distances between the stations could be uniquely deter- 
mined. This discovery gave rise to the idea of fixing one or more baselines 
between the stations in the planar four -station system. Analysis showed that 
this four-station system did not degenerate mathematically once two of the base- 
lines were fixed. Hence, a workable system could be obtained by aligning sta- 
tions in straight lines and using the mathematically stable straight line solutions 
to determine the baselines which then would be used to yield the numerically 
stable four station solution. Such a system is actually possible for short 
baselines. 

More analysis on the geometry peculiar to this problem showed that the 
addition of two more stations was equivalent to fixing the ground baselines 
between two sets of stations. Therefore, as was expected and further theoretical 
and numerical analysis indicated, a six-station configuration permitted simul- 
taneous determination of the station coordinates and the position/velocity of the 
satellite if at least four simultaneous fixes or strikes are taken by the station 
net, even if the stations are nearly coplanar. Subsequent analysis has shown 
that accuracies in station coordinates to a centimeter level can be obtained. 
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Lately, it has been demonstrated that the systems described herein can 
also function if only range differences, instead of range only information is 
available. 

Further analysis will no doubt uncover more information relative to the 
peculiarities of such geometric ranging schemes. Nevertheless, the feasibility 
of the geometric method for station location determination and ephemeris com- 
putation has been clearly established. 
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SECTION III 


INTRODUCTION AND COMPARISON TO OTHER SYSTEMS 

^ In this section other techniques will be discussed which might be used to 
monitor tectonic motion and implement earthquake hazard assessment. The 
other techniques to be considered are the Geodolite, Tiltmeters, Gravimeters, 
very long baseline interferometry (VLBI), and the Goddard Space Flight Center 
laser ranging system. The salient features of each system will be described 
first, and then a comparison of the described systems with the 3-D Multilatera- 
tion technique proposed in this report will be presented. 

3. 1 THE GEODOLITE 

The Geodolite is a laser- ranging device which is capable of measuring the 
distance between two points on the Earth's surface, but is not capable of resolv- 
ing this distance into 3-dimensional components. The device is manufactured 
by Spectra Physics (Mountain View, California) and is currently being used by 
the U.S. Geological Survey. The Geodolite transmits an amplitude -modulated 
laser beam to a ground-based reflector; the modulation of the reflected beam is 
then compared with the modulation of the transmitted beam in order to evaluate 
the range between transmitter and reflector. The U.S. Geological Survey cal- 
culates the distances between benchmarks (reference points on the ground) by 
grouping the benchmarks into sets of 7, and then using a Geodolite to measure 
each of the 15 distances between different pairs of the 7 benchmarks. Since the 
measurements are partially redundant some statistical averaging is possible, 
and inter station distances up to 30 km can be evaluated with a standard deviation 
of error which is less than 1 cm. Clearly, the Geodolite provides excellent 
accuracy over its severely constrained range. However, in order to achieve 
the quoted accuracy, it is necessary to correct the crude data with temperature 
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and humidity measurements obtained by flying a helicopter along each path 
that is measured by the Geodolite. 

Although it is conceivable to combine a large number of Geodolite mea- 
surements in order to measure distances up to 1000 km, there will be some 
serious problems. For example, if a distance of 900 km is partitioned into 
30 segments of 30 km each, and each segment is measured by using the Geodolite 
and the U.S. Geologic Survey technique then: 

• It will be necessary to measure 450 (30 x 15) separate distances. 

• The overall standard deviation of error will be at least: 

\, r 30 = 5.5 cm 

• To obtain the complete set of measurements will require a great deal 
of time and rapid crustal deformations (which may give premonition 
of earthquakes) will be lost. 

Furthermore: 

• The Geodolite is practical only when measuring distances between 
points which are elevated with respect to the intervening terrain. 
(Otherwise the line-of- sight distance will be reduced to a few 
kilometers, and variations in ground temperature will cause 
anomalous propagation of the light beam). It is extremely unlikely 
that the 450 required measurement paths can be chosen so that 
each path connects points which are elevated with respect to the 
intervening terrain. 


*This assumes that the errors in measuring each of the 30 segments are uncor- 
related. If the errors are correlated, the overall standard deviation will be 
greater than 5. 5 cm. 
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• Since the Geodolite can only measure lateral distance, vertical 

movement of the benchmarks will not be easily discernible. How- 
ever, many important strain fields (such as the field related to the 
1971 San Fernando earthquake in the Los Angeles area) have verti- 
cal components. 

3. 2 TILTMETERS AND GRAVIMETERS 

Tiltmeters and gravimeters could conceivably be used in conjunction with 
a Geodolite so as to provide independent measurements of vertical motion. 

Both tilmeters and gravimeters have the capability of detecting vertical ground 
motions within the order of centimeters. Tiltmeters show only a change in 
slope, however, and not the amplitude of a vertical ground excursion at any 
one place. Tilt measurements are thus not directly comparable with the 
results anticipated from laser retroreflection. 

Gravimeters, on the other hand, can detect the amplitude of ground 
motion because of the change in gravity which results from a change in distance 
from the center of the Earth. The change in gravity is equal to the magnitude 
of the "Bouguer correction" for gravity analysis (see Glossary). It is the 
change in gravity due to a change in elevation when the ground surface moves 
along with the instrument. 

The accurate value of the Bouguer correction at any specific location 
depends partly upon the local rock density, but it is about 2 microgals. ^ 

The modern LaCoste -Romberg gravimeters have a sensitivity capsule of 
reading to one microgal, accuracies of 3 pgal with repeated exceedingly care- 
ful observations, and 10 pgal with ordinary good work. The wideband 


"'One pgal is an acceleration of 10 


-6 


2 -9 

cm /sec and about equal to 1 0 g. 
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accelerometer described by Block and Moore [3. l] also has a sensitivity and 
accuracy capability of measuring elevation changes of about one centimeter. 

Both instruments have a considerable drift rate, however, which is the 
equivalent of several centimeters per week. Neither instrument could be left 
by itself to monitor elevation changes at individual stations because the drift 
rate would make the results uninterpretable at worst and ambiguous at best. 

In order to account for instrumental drift rates, a gravimeter program 
would involve frequent transportation of the reading gravimeter to a base station 
for calibration, or the continuous transportation of one gravimeter from a base 
station to the reading stations. Either case would require constant work by a 
ground staff. 

In summary, neither tiltmeters nor gravimeters appear practical for 
measuring vertical motion to the required accuracy. 

3. 3 THE GODDARD LASER RANGING SYSTEM* 

Goddard Space Flight Center has designed and built a laser ranging sys- 
tem which employs two ground stations and which calculates the geocentric 
coordinates of one station given the geocentric coordinates of the other station. 
Each station transmits a laser pulse to a retror eflecting satellite, receives the 
reflected pulse, and computes station-to-satellite range from the measured 
value of time between transmission and reception of the laser pulse. The range 
data from each station is then combined with 1) an extremely detailed model 
of the forces acting on the satellite, and 2) the assumed geocentric coordinates 
of one station, in order to yield the geocentric coordinates of the other stations. 


' The information presented in this section is based on private communication 
from Dr. David Smith of GSFC as related to R. Jaffe. 
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Goddard has tested its system by ranging to satellites at altitudes of 
1000 km from stations which are separated by 400 km. The station-to- satellite 
ranges were measured with a standard deviation (about the measured mean) of 
about 3 5 cm. Inter station distance was evaluated with a standard deviation of 
approximately 25 cm, and the 3 -dimensional components of the inter station dis- 
tance were evaluated with a standard deviation of approximately 50 cm. 
Goddard's error analysis indicates that inter station distances were calculated 
with an accuracy of about 2 meters, and that the 3 -dimensional components of 
interstation distances were calculated with an accuracy no better (but not 
worse) than 2 meters. 

Goddard expects to be able to improve the accuracy of its system by 
utilizing more advanced laser -ranging components, and more sophisticated 
predictions of the force fields on the satellite. It should be noted that Goddard 
already employs an extremely sophisticated trajectory program which includes 
the effects of polar motion, Earth tides, and gravity anomalies. Goddard esti- 
mates that its present orbit determination model is accurate to a few meters 
over a 2-3 day period for a 1000 km altitude satellite. However, the present 
accuracy of the Goddard technique must be improved by an order of magnitude 
before the system can be used in the earthquake -hazards applications described 
in this report. 
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3.4 VERY LONG BASELINE INTERFEROMETRY 


This technique [3. 2] has its origins in optical and radio astronomy. The 
technique operates by receiving random radio signals generated by objects 
external to our own galaxy, typically, quasars and other galaxies. By measur- 
ing the difference in time of arrival of identical random signal patterns at a 
pair of stations, the distance between the stations can be determined. 

An error frequently made in projecting the usefulness of this system to 

geophysical goals is to confuse the angular resolving power of the interferometer 

with an equivalent linear distance on the Earth's surface. Since the resolving 

power of an interferometer is given by the ratio of the received wavelength to 

- 8 

the distance between the receiving stations, resolutions of 10 arc seconds 
are possible. It is not correct to conclude that such resolution alone provides 
Earth distances to be measured with accuracies of 1 centimeter. 

The ability to invert very long baseline interferometry (VLBI) observations 
into accurate geophysical measurements is limited mainly by Earth atmospheric 
uncertainties. However, dual frequency VLBI holds significant promise for 
ionospheric calibrations and water vapor radiometers appear promising for cali- 
brating the wet component of the troposphere. Therefore, VLBI performance 
at the 10 cm baseline accuracy level is being forecast within a few years and 
5 cm or better for the long term. 

VLBI possesses the following advantages: 

• Capacity to determine three-dimensional station coordinates in a 
geocentric coordinate system accurate to 10 cm. 

• Capacity to operate under virtually all weather conditions. 

• Measurements are made relative to a frame of virtually time 
invariant extragalactic radio sources. 
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Portable stations are feasible for remote station utilization. There 


is no VLBI requirement to build or launch any spacecraft or to 
transmit from Earth any radio or light signals. VLBI is entirely 

passive, using only natural radio signals. 

Some disadvantages are: 

• Water vapor in the atmosphere prevents precise system calibration. 

• Software costs are relatively high. 

• Some modeling errors are present. 

It is interesting to speculate how the present multilateration technique and 
VLBI might be combined (use of range difference data from a satellite) to over- 
come some of the inherent disadvantages of each method. 

3. 5 COMPARISON OF THE 3-D MULTILATERATION TECHNIQUE WITH 

THE DESCRIBED SYSTEMS 

The Geodolite appears to be suitable for many of the applications envi- 
sioned for the 3-station collinear system. Both the Geodolite and the 3-station 
system can measure short distances (but not the 3 -dimensional components of 
these distances) with high accuracy and at relatively low cost. The Geodolite has 
the advantage of requiring considerably less capital investment; however, the 
Geodolite can measure only line -of- sight distances, and can achieve high 
accuracy only when its raw measurements are corrected with temperature and 
humidity data obtained by overflying the baseline with a helicopter. Neither 
the Geodolite nor the 3-station system are applicable to the most difficult mea- 
surement problem in earthquake hazard assessment, viz. , accurate evaluation 
of the 3 -dimensional components of strain between widely separated stations. 

Tiltmeters and gravimeters initially seemed to offer some promise of 
augmenting the distance measurement capability of the Geodolite with inde- 
pendent measurements of vertical motion. However, more detailed investigation 


JPL Technical Memorandum 33-605 


3-7 



revealed that tiltmeters were incapable of directly measuring either distance 
or motion, while gravimeters were incapable of yielding the requisite accuracy 
without frequent recalibration. 

The Goddard laser ranging system shares many features in common with 
the 4- and 6-station systems. All systems employ laser ranging, and all sys- 
tems are capable of measuring the 3 -dimensional components of distance 
between widely separated stations. 

The Goddard system has the advantage of requiring fewer stations, and 
hence requiring a smaller capital investment for a complete operating system. 
However, Goddard needs approximately 2 months of range measurements to 
yield accurate output data, whereas JPL systems can yield accurate output data 
in only hours within reception of the ranging measurements. Consequently, if 
many locations are to be mapped in a limited period of time, only one of the 
JPL systems would be needed (since data could be obtained rapidly and the 
stations then transported to new locations), whereas a number of Goddard sys- 
tems would be needed. The systems' rapidity of operation is also significant 
in applications such as detection of sudden ground deformations which often 
give a premonitory indication of earthquakes, and monitoring of seismic acti- 
vity during the period immediately following a major earthquake. 

Although cost, complexity and time required for observation are all 
important bases for comparison among systems, the ultimate criterion for 
system evaluation, at least with regard to earthquake hazard applications, is 
accuracy. And it is in terms of accuracy that the Goddard system and the 
3-D Multilateration technique manifest their most significant differences. The 
accuracy of the Goddard system is fundamentally dependent upon the accuracy 
of the satellite force-field model; unless the accuracy of this model is 
improved 10-fold or more, it is extremely unlikely that Goddard will be able 
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to achieve station location accuracies of 10 cm or better, regardless of any 
improvements made in their laser ranging equipment. In contrast, the 3-D 
Multilateration technique is essentially independent of Earth force fields, and 
is constrained only by the accuracy of the laser ranging equipment. Since 
excellent laser ranging accuracies can be achieved using commercially available, 
hardware, the 3-D Multilateration technique will be able to achieve station 
location accuracies of 1 cm or better. It seems highly improbable that the 
Goddard system will be able to match this performance in the foreseeable future. 

On the other hand VLBI techniques which use the principle of multilatera- 
tion, wherein the raw data is the time delay between reception of a signal at two 
stations from a satellite (which is subsequently transformed into the slant range 
difference between the satellite and each of the two stations) seem to offer a 
second avenue of approach. Specifically, the geometric principles invoked in 
this report can be transformed such that range differences instead of range 
itself become the raw data. Under these transformations it is possible that 
VLBI might be able to reach to within a few centimeters of the accuracy of the 
laser ranging system. 
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SECTION IV 


COORDINATE SYSTEMS FOR RANGING EXPERIMENTS 

In this section, the inertial coordinate system of dynamical astronomy 
is outlined in order to provide the reader with a means of comparing the coor- 
dinate system to be introduced in this study with the standard system. Deter- 
mination of station coordinates relative to an inertial system briefly is dis- 
cussed next. Geometric coordinate systems peculiar to laser ranging are dis- 
cussed and a comparison to the inertial system is made. Finally, the funda- 
mental coordinate system utilized in this study is defined. It is concluded that 
the geometric coordinate system is a self-consistent system capable of yield- 
ing station locations with high accuracy relative to a plane defined by three of 
the stations. The accuracy with which the true equator can be determined 
relative to the geometric equator is estimated, and explicit formulas for the 
transformation of coordinates are discussed, 

4. 1 INERTIAL SYSTEMS 

Newtonian laws are always referred to a conveniently adopted inertial 
coordinate frame. The system commonly utilized in dynamical astronomy is 
displayed in Figure 4-1. 

In this coordinate system the principle axis, x, is taken to pass 
through the true equator and equinox of date. This means that the position of 
the pole of the Earth has been modeled as accurately as is presently permitted, 
i.e. , all the precession and nutation effects caused by the perturbative action 
of the Sun and the Moon on the equatorial bulge of the Earth have been incorpo- 
rated into the coordinate system model. The y axis is in the true equator 
(where the equator actually is at a given instant) and the z axis is taken to 
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Figure 4-1. 


Inertial Coordinate System 
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complete the right handed coordinate system. Thus, the z axis is directed 
along the best possible estimate of the Earth's angular momentum vector. 

In a broad sense, if the wandering of the pole caused by the lack of 
alignment of the Earth's body fixed system and the Earth's angular momentum 
vector is neglected, the linkage between a geographically fixed system and the 
best estimate of an inertial system is the true sidereal time . 

On Earth, time is usually reckoned in a mean sense, e. g. , as mean 
sidereal time. To determine the actual or true position, the true sidereal 
time as opposed to the mean sidereal time must be determined. To determine 
the true sidereal time, the equation of the equinoxes [4. l], [4. 2], must be 
evaluated. This equation is defined as 

E = 64* cos + 6t), 


where 6^ is the nutation in longitude, 6«= is the nutation in obliquity, and is 
the mean obliquity of the ecliptic. From observational astronomy the funda- 
mental constants of these expressions can be obtained, and analytic expressions 
with time as the independent variable can be derived based upon theory [4. 3], 

[4. 4j. Once E is evaluated, the true sidereal time, 0^, is given as 


0 g - 0 T - 0 M + E > 

where 9^, the mean sidereal time, is given as a polynomial in time [4. 1]. 

The local sidereal time at any observational station is then taken to 


be 


9 - 9rp + X.£ > 

where X._ is the east longitude of the station. The rectangular coordinates of 

hj 

the station from the dynamical origin are then defined in Appendix J. 
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DETERMINATION OF STATION COORDINATES 


4. 2 

If the study objective is to use accurately measured data, e. g. , laser 
range measurements, in order to determine the locations of ground stations, a 
few fundamental points must be kept in mind. 

First, laser data is taken in a coordinate system with respect to the 
body-fixed axes of the Earth at any given instant. Second, the satellite position 
will be computed in the true inertial coordinate system defined previously. 

Hence, if the satellite ephemeris is used to determine a station coordinate, say, 
the latitude, that coordinate will be determined relative to the inertial coordinate 
system. 

Therefore to compare the new estimates of station coordinates with the 
old, the rectangular coordinates of the stations should be transformed back to 
the body-fixed axes to which geodetic coordinates are referred. The trans- 
formations are standard, [4. 1], [4. 5| . 

However, since the Conventional International Origin of the Earth's 
body-fixed coordinate system was actually chosen arbitrarily [4. 6], it follows 
that, if a new method for estimating station coordinates of much higher accu- 
racy becomes available, it might be desirable to adopt a new equatorial position 
and reference all new station data to such a new benchmark. If the new method 
should be much better than the old method, a much more consistent set of 
station coordinate data would be available. 

4. 3 GENERALIZED GEOMETRIC COORDINATE SYSTEM 

It will be demonstrated in this study that, if simultaneous range 
measurements are taken from a suitable number of ground stations to an orbiting 
satellite, then the rectangular coordinates of the stations can be determined to 
extremely high accuracy in a suitable coordinate system without knowledge of 
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the trajectory. Without going into further details let the following coordinate 
system be examined (Figure 4-2). 

In Figure 4. 2, the geometric equator is defined by physically position- 
ing two stations to be on the best possible estimate of the equator on a given 
adopted date. By then adopting the coordinates of these stations by definition to 
be [X 2 , 0, o], [X 3 , Y 3 , o] , it follows that the geometric equator on this date 

has been defined relative to some origin [0, 0, Oj. Hence the X axis passes 
through X 2 , the Y axis lies in the plane defined by stations 2 and 3, and Z is 
perpendicular to the plane defined by the X-Y plane. 



Figure 4. 2. Geometric Coordinate System 
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Since current estimates of station position with respect to the equator 
can only be accomplished to within a few meters, it follows that the dynamical 
coordinate system and the geometric system can only be brought into agreement 
by an error source of this order of magnitude. 

It will be seen that, via laser ranging, the coordinates jjX^, X^, Y^> 

X^, Y^, Z^, X^, Y^, can be computed within the geometric framework 

two orders of magnitude more accurately than in the dynamical coordinate sys- 
tem. Furthermore, once the coordinates of the stations are obtained, it 
follows that the satellite coordinates [x, y, z]at that instant can be derived. 

The negative of these coordinates is therefore the location of the origin of the 
geometric system as viewed by the satellite. 

In summary, the geometric system yields station coordinates in a 
relative framework, in which the accuracy of relative distance measurement is 
only limited by laser equipment accuracy. 

4. 4 ADOPTED GEOMETRIC TOPOCENTRIC COORDINATE SYSTEM 

The previous discussion indicated that to make the generalized geo- 
metric system agree with the standard astronomical system, the longitude of 
station Z with respect to the Greenwich meridian would have to be precisely 
known. Furthermore, stations Z and 3 would have to lie precisely on the equa- 
torial plane. Since exact alignments of this sort are not possible, it makes 
sense to adopt some well known surface station as the origin of a topocentric 
system to be used in laser ranging studies. This origin will of course have 
errors in latitude, longitude and elevation relative to the inertial astronomical 
system of a few meters, i. e. , a bias-like error in the coordinates with respect 

Ox 

to the adopted origin of the dynamical center. The adopted topocentric system 


The origin of the dynamical system cannot be defined to any greater accuracy 
than the origin of the generalized geometric system. 
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Y 

Figure 4-3. Geometric Topocentric Coordinate System. 


is defined in Figure 4-3. As can be seen, the origin is placed at some 
geodetically known reference station, R, i. e. , Station 1 whose geodetic 
coordinates are k, 4> , H. The second station defines the X axis; the third sta- 
tion defines the XY plane; the Z axis is normal to the XY plane and is directed 
away from the Earth, and the Y axis completes a right handed system. This 
coordinate system will be used throughout the present report and is referred to 
as the adopted coordinate system. 

4. 5 STATION COORDINATE RELATIONSHIPS 

For purposes of systems analysis, specifically for actually computing 
sets of realistic station coordinates in the adopted geometric coordinate system 
from tabulations of geodetic station coordinates which are presently known, a 
transformation of coordinates will be required. This transformation is explicitly 
developed in Appendix J. 
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MULTI LATER ATION: 

A METHOD OF MEASURING INTERCONTINENTAL 
PLATE MOVEMENTS 
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SECTION V 


EXPECTED SYSTEM ACCURACIES FOR REAL WORLD CONFIGURATIONS 

£ 

JpZf In this Section, numerical results are presented pertaining to the 
system accuracy achievable by use of the geometric principles discussed in 
Section I. 

To illustrate the versatility of the primary recommendation of this 
report, i. e. , use of the 6-station system, four hypothetical experimental con- 
figurations using six stations are evaluated, namely: 

• an intermediate distance network along the Southern California 
portion of the San Andreas Fault, 

• a long-baseline network covering the continental United States, 

• intercontinental networks of stations, 

• a very short baseline network within the Los Angeles metropolitan 
area. 

As will be seen, each of these configurations yields accurate solutions for the 
station coordinates which , subject to the recommendations of Paragraph 6.4, will 
ultimately approach accuracy levels of one centimeter. 

The accuracies of the four- station, intercontinental baseline system 
and the three- station, short baseline system also are discussed briefly. 

It should be noted that the choice of stations and satellite (or airplane) 
trajectory points to be used in the following applications is rather arbitrary, 
although in a few cases selection has been guided by preliminary geophysical 
considerations. No attempt has been made to achieve optimal configurations; 
however, care has been taken to avoid certain degenerate configurations (see 
Section VII). It is believed that with suitable optimization techniques (a task 
which is proposed in Section X for future study), and with the use of more 
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trajectory points than those employed in these calculations, even more dramatic 


accuracies than presently reported will be achieved. In all tabulations to fol- 
low, the average errors, € , in each station coordinate and the standard devia- 
tion, or, have been computed on the basis of simulation results for 10 different 
sets of random range errors, where each set consists of 100 points drawn from 
a Gaussian distribution with zero mean and standard deviation of 1 cm. * Such a 
procedure may be interpreted in the following manner: assume that 100 strikes 
of the satellite(s) per pass are possible; then the average errors quoted are 
what would be expected to be obtained by averaging the results over 10 identical 
passes; furthermore the standard deviations are the standard errors of the cal- 
culated station locations from one typical pass of data. 

The procedures for calculating the various quantities used in the follow- 
ing subsections are outlined in Section VII, and detailed computational algo- 
rithms are given in Appendix A and Appendix C. 

5. 1 SIX-STATION CONFIGURATION: A SYSTEM FOR ALL BASELINES 

The use of six stations is the primary recommendation of this report 
(see Section X). A workable six- station system conceivably can be set up in 
any global or local region of geodetic interest. In contrast, the four- station 
and the collinear three- station configurations (as discussed before) will yield 
accurate station locations only for extremely wide and short baselines, respec- 
tively. Examples of the kind of accuracies that can be obtained in a six-station 
configuration are given in the following subsections, 

5. 1. 1 Local Experimental Network 

Figure 5-1 shows the locations of six hypothetical stations placed along 
a region of great interest to geophysicists and seismologists, i.e., the Southern 

^Precise definitions of the quantities € and cr are given in Paragraph 7. 3. 

System Bias errors are discussed in Appendix K. 
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California portion of the San Andreas Fault. The coordinates of these stations* 
in the adopted coordinate system (see Paragraph 4.4) with San Simeon as the 
origin are displayed in Table 5-1 and Figure 5-1. 

Table 5 - 1 . 


Adopted Coordinates in Kilometers of the 
San Andreas Fault Stations** 


Station 

X 

Y 

Z 

1 . 

San Simeon 

0 

0 

0 

2. 

San Diego 

495 

0 

0 

3. 

Isabella 

186 

160 

0 

4. 

Santa Rosa Island 

172 

89 

+ 5. 1 

5. 

Millerton 

6. 2 

200 

- 6. 1 

6. 

Blythe 

599 

235 

-12. 9 


'■'^hese coordinates have been obtained from a simple computer program 
accepting as input the geographic longitudes and latitudes of the stations. 

A spherical Earth was assumed. 

**The Z coordinates are apparently large due to the adopted coordinate system, 
wherein stations 1 and 2 define a chord cutting the globe; this has the effect 
of placing, e. g. , station 6 at a sizeable distance from the X-Y plane. Note 
that the scale in the sketch has been exaggerated. 
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Two circular orbit satellite passes, whose projections on the XY plane 
of the adopted coordinate system are indicated in Figure 5-1, were used in the 
simulation. The heights of the two satellites are 500 km and 750 km for orbits 
1 and 2, respectively. Fifty points from each orbital arc (a total of 100 points) 
were generated with the X-coordinate s ranging from San Simeon to Blythe. The 
accuracy simulation procedure as outlined in Paragraph 7. 3 was then performed. 
Table 5-2 summarizes the results of the simulation for this system configuration. 

Table 5-2. 

The Average Errors, e , and the Standard Deviations, cr , 
of the Coordinates of the San Andreas Fault Stations 
in Centimeters Assuming 1 cm 
Normal Random Errors in Ranging 


Station 

€ 

X 

(T 

X 

6 Y 

cr 

Y 

£ z 

°Z 

2. 

San Diego 

0. 15 

2. 06 





3. 

Isabella 

-0. 15 

1. 11 

0. 84 

3. 56 



4. 

Santa Rosa Is. 

-0. 20 

1. 31 

-0. 05 

1. 76 

0. 10 

0. 99 

5. 

Millerton 

0. 09 

1. 14 

0. 91 

4. 16 

0. 03 

0. 71 

6. 

Blythe 

0. 26 

2. 64 

1.41 

4. 74 

0. 24 

1. 18 


5.1.2 Continental Network 

A hypothetical six- station configuration spanning the continental United 
States is shown in Figure 5-2, with the adopted station coordinates displayed in 
Table 5-3. The assumed satellite trajectory data is also shown on the figure. 

The numerical simulation resulting for this configuration is summarized 
in Table 5-4. 
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Figure 5-2. The Continental Station Network 


Table 5-3. 


Adopted Coordinates in Kilometers 
of the Continental Stations 



X 

Y 

Z 

1 . 

Goldstone, Cal. 

0 

0 

0 

2. 

Cape Kennedy, Fla. 

3740 

0 

0 

3. 

Spokane, Wash. 

120 

1350 

0 

4. 

Denver, Col. 

1060 

510 

250 

5. 

Houston, Tex. 

2190 

-82 

270 

6. 

Washington, D. C. 

3310 

1400 

75 


Table 5-4. 

The Average, € , and the Standard Deviation, cr , 
of Station Coordinate Errors in Centimeters 
for the Continental Stations, Assuming 1 cm 
Normal Random Error in Ranging 


Station 

€ x 

T x 

€ Y 

r Y 

e Z 

r Z 

2. 

Cape Kennedy 

-0. 85 

3.02 





3. 

Spokane 

0. 21 

1. 70 

0. 54 

3.81 



4. 

Denver 

-0. 37 

1. 32 

-0.09 

1. 22 

Q. 25 

0.52 

5. 

Houston 

-0. 39 

1. 59 

-0.28 

1.02 

Q. 32 

1.00 

6. 

Washington 

-0. 25 

3. 84 

0.39 

3. 90 

0. 32 

1. 35 
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5. 1. 3 Intercontinental Networks 


The determination of the coordinates of the stations located on different 
continents to centimeter accuracy is of great interest to geologists, tectonic 
physicists and astronomers. * 

It should be remarked that less-than-10 cm accuracy station location 
can be achieved with only four intercontinental stations, provided that the four 
stations are favorably located. With six stations, however, computational expe- 
rience indicates that obtaining accurate solutions seems to be the rule rather 
than the exception. Further analysis will explicitly define this rule. 

As an illustration, consider Figure 5-3 where a hypothetical configura- 
tion is adopted at the following cities: New York (U. S. A. ), London (England), 
Mexico City (Mexico), Caracas (Venezuela), Freetown, Sierra Leone (Africa), 
and Reykjavik (Iceland). A total of 100 trajectory points from two satellites at 
altitudes 6, 000 and 10, 000 kilometers, respectively, were employed in the cal- 
culation. The accuracy corresponding to this configuration is displayed in 
Table 5-5. 

As a second illustration, the following choice of sites were selected: 
Kaneohe, Hawaii (U.S.A.), Tarara (Peru), Calgary (Canada), Mexico City 
(Mexico), San Bernardino (U.S.A.), and Clipperton I. Two circular-orbit satel- 
lites of 6, 000 and 10, 000 km altitudes also were employed. Table 5-6 sum- 
marizes the results of the accuracy simulation. 


*With this kind of accuracy, the hypothesis of continental drift will be subject 
to direct scrutiny. 
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Figure 5-3. An Intercontinental Network 
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Table 5-5 


The Average, e, and the Standard Deviation, tr , 
in Centimeters of the Errors in the 
Station Coordinates of the 
Intercontinental Configuration 
(Origin Assumed to be 
at Mexico City) 


Station 

€ x 


n 

cr y 

€ Z 


London 

-0. 94 

4.06 





New York 

-0. 82 

2. 39 

-0.83 

2. 07 



F reetown 

1 

o 

• 

oo 

4. 10 

0.34 

2. 83 

0.45 

1.66 

Reykjavik 

-0. 77 

3.60 

o 

m 

• 

o 

1 

0. 89 

0.02 

0. 41 

Caracas 

-0. 17 

1.64 

-0. 26 

2. 35 

0. 18 

0. 94 


Table 5-6. 

The Average, e, and the Standard Deviation, cr , 
in Centimeters of the Errors in the 
Station Coordinates of the Second 
Intercontinental Configuration 
(Origin Assumed to be 
at Kaneohe) 


Station 

€ X 

°x 

£ y 

r Y 

e z 

^Z 

Talara 

1. 06 

3. 84 





Calgary 

2. 12 

4. 41 

-.37 

4. 16 



Mexico City 

1. 09 

3. 20 

-.60 

1. 16 

. 19 

. 34 

San Bernardino 

1.41 

2. 63 

-.40 

2. 08 

-.00 

. 54 

Clipperton I. 

0. 4.9 

2. 25 

-.60 

.93 

-.02 

. 73 
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5.1.4 Airplane Experiment 


In extreme contrast to the previous experiment, consider a small scale 
experiment that can be performed within the Los Angeles basin with the stations 
depicted in Figure 5-4. An airplane was used in lieu of satellites. Because 
of the short baselines, the stations are located almost on a plane with vanishing 
Z -components ; however, this does not affect the accuracy of the six-station sys- 
tem. The usual assumption about the ranging accuracy (1 cm) has been imposed, 
Two airplane trajectories, as indicated in Figure 5-4, were used in the simula- 
tion. The resulting accuracies of various station coordinates are shown in 
Table 5-7. 


Table 5-7. 

Average, « , and Standard Deviation, cr , in Station Locations 
for the Airplane System in Los Angeles 
(Origin Assumed to be at San Fernando) 


Station 

€ X 



Ty 

£ Z 

17 Z 

Pasadena 

-1.28 

1. 57 





Santa Monica 

4.05 

3.87 

0. 39 

2.24 



Burbank 

-0. 13 

0.48 

0.38 

0. 32 

0.49 

Q. 91 

Mt. Wilson 

-0. 70 

1.67 

0.39 

0.47 

0. 38 

1.25 

V an Nuys 

1. 37 

1. 33 

0.30 

0.75 

-0. 20 

0.77 
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FOUR STATION CONFIGURATION: INTERCONTINENTAL BASELINES 


5. 2 

As mentioned in the summary (Section I) and will be discussed in 
Section VII, a system of four stations, each ranging to at least six trajectory 
points in principle is a workable configuration provided that the four stations do 
not lie on a common plane. (Certain orbital degeneracies must also be ruled 
out; see Appendix B.) However, owing to the Earth's flatness over small dis- 
tances, the four- station system can be utilized only if baselines are of intercon- 
tinental scale and if one of the stations is separated by a sizable amount from 
the plane containing the other three stations. Even if these two conditions are 
satisfied, there is no guarantee that the system will yield coordinates accurate 
to the centimeter level. However, there do exist four- station configurations 
that give acceptable accuracies. An example of this is the following choice of 
sites: 

Station 1: Mexico City 
Station 2: London 
Station 3: New York 
Station 4: F reetown 

Note that the above configuration is a subset of the six- station intercontinental 
network discussed in Paragraph 5. 1. 3. The same satellite trajectories as indi- 
cated in Figure 5-3 are also used in this analysis. The results of the accuracy 

analysis for this configuration are shown in Table 5-8. 

Table 5-8. 

Average, e, and Standard Deviation, tr , in Centimeters of Station Coordinate 
Errors for the Four-Station Intercontinental System 
(Origin Assumed to be at Mexico City) 


Station 

e X 

7 x 

€ Y 

Y 

£ z 


London 

-3.93 

7. 78 





New York 

-1.77 

2.22 

0. 51 

5.91 



F reetown 

-4.02 

13.49 

1.83 

2.25 

1.28 

1.63 
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A 4- station simulation was also carried out with stations situated at 
Kaneohe, Tarara, Calgary and Mexico City. Comparable results to the 4- 
station intercontinental simulation were obtained but are omitted here. 

As can be seen by comparing Table 5-8 with Table 5-5, the accuracy 
of the four-station system is inferior to the corresponding six- station system. 

5.3 COLLINEAR THREE -STATION CONFIGURATION: VERY SHORT 

BASELINES 

If three stations are aligned in a straight line, then there exists a 
simple technique for surveying the baselines between the stations using only two 
simultaneous range measurements to an aircraft. Obviously, such a technique 
can be applied only to short baselines, since it is impossible to align widely 
separated stations in a straight line due to the curvature of Earth's surface. A 
discussion of the collinear three- station technique is presented in Paragraph 7. 5 
and in Appendix C. 

The accuracy of the collinear three- station system depends on how 
well the stations can be aligned in a straight line and how the trajectory points 
are actually chosen. Let be the amount of misalignment of station 3 with 
respect to the straight line (X-axis) joining stations 1 and 2, as illustrated in 
Figure 5-5. 


Y 



Figure 5-5. Misalignment of Station 3 in the 
Collinear Three-Station System 
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It can be shown that the trajectory points can be chosen in such a way 
that the induced errors in one of the station coordinates (X^ or X^) due to the 
misalignment error can be minimized. However, the corresponding errors 
in the other station coordinate (X^ or X^) will remain large.* This implies 
that each of the two baseline X^ and X^ should be determined separately with 
suitably- chosen trajectory points. 

For illustration purposes, the following nominal baselines are adopted: 
X^ = 40 km 
X^ = 60 km 

Let the retroreflector-bearing vehicle fly a planar trajectory given by the 
equation 

y - 0. 5(x - X 3 ) 

For various degrees of misalignment Y^, the induced errors in X^ and X^ are 
computed, with the usual assumption of normally distributed ranging errors 
having zero mean and 1 cm standard deviation. The use of 100 trajectory points 
is also assumed. Results of this calculation are summarized in Table 5-8. 

As can be seen from this Table, the errors in X^, with misalignment 
error Y^ as large as 100 meters, remain well within the centimeter level even 
though errors in X^ become progressively larger as Y^ increases. 


*See Appendix C for more details of the error minimizing technique. 
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Table 5-9. 


Accuracy of the 3-Station System as a Function of the 
Misalignment from the X-axis. Trajectory points 
have been chosen to minimize the induced error in 2C ? . 
Ranging errors are Normally Distributed with Zero 
Mean and 1 cm Standard Deviation. The Average Values 
for 10 Calculations each using 100 Trajectory Points are 
Quoted. All quantities are in Centimeters. 


Misalignment Y^ 

Error in 

Error in X^ 

1 

0. 287 

0. 683 

10 

0. 101 

5. 10 

100 

0. 231 

50. 2 

1, 000 

0. 020 

500. 

10, 000 

0. 350 

5000. 


Similarly, the errors in can be minimized by flying the retroreflector on a 
planar course satisfying, e. g. , 
y = - 0. 5( x - X 2 ). 

In such a case, the errors in X^ would remain small while the X^ error increases 
with increasing misalignment. 

5.4 ULTIMATE SYSTEM ACCURACY 

It must be understood that the results quoted in this section are pre- 
liminary and can be improved. In the historical development of this project 
(see Section II) it was pointed out that the six-station configuration was a recently 
developed concept. Due to this fact and lack of time, it has not been possible to 
determine the ultimate system accuracy. This accuracy depends on two factors, 
namely, the manner in which orbital orientations are adopted and upon station 
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location. The first factor is flexible and upon proper orbital design, accuracy 
improvements will be noted. The second factor is less flexible due to geographic 
constraints, but certain optimization can still be performed. These two factors 
coupled with: 

• the use of more data, since only 100 measurements were utilized 
in most of the simulations, and 

• the use of additional stations, 

will result in station coordinate determination whose accuracy is only limited by 
the accuracy of the laser measurements, i. e. , the mapping of errors from the 
measurements to the station coordinates will be better than one-to-one. There- 
fore, since laser measurements accurate to approximately 1 cm can be obtained 
presently (see Section X), the following conclusion can be drawn: A laser 
ranging system accurate to the sub- centimeter level can be constructed within 
the next five years. 
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SECTION VI 


GEOLOGICAL APPLICATIONS OF MULTILATERATION: 

STATION LOCATION SELECTION 

C 

A major geophysical objective of accurately locating widely spaced points 
as a function of time is to monitor the strain buildup in the Earth's crust. Then, 
comparing the strain buildup with the strain release resulting from earthquakes 
and from nonseismic slippage or creep, a model of elastic deformation and its 
variation from place to place can be constructed. In more detail there are three 
essential components in the determination of earthquake hazard to be anticipated 
in the major earthquake zones of the world. These three features include: 

(1) The strain field (see Glossary) of the crust of the Earth, which describes the 
state of deformation at any given time and how this field changes with time; 

(2) the stress field (see Glossary) in the Earth's crust and how this field changes 
with time, and (3) the critical states for rupture and tearing along earthquake 
faults . 

In this section the critical components required for earthquake estimation 
will be identified. This will be followed by a brief discussion of plate tectonic 
theory. These discussions will naturally lead to a number of pertinent ques- 
tions, Some of these questions, specifically those relative to the strain field 
in the Earth's crust will no doubt be answered by means of the proposed multi- 
lateration techniques. 

The analysis in this Section is intended to be of a general nature and 
depicts the process of earthquake prediction as a complex process invoking 
many factors. For another viewpoint to this problem the reader is directed to 
Appendix L. 
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6. 1 CRITICAL COMPONENTS OF EARTHQUAKE ESTIMATION 


The problems of earthquake prediction can be most simply discussed by 
reference to Figure 6-1. Imagine that a sample of rock or other material is 
being loaded under a constant load stress, cr . Let it be assumed for the moment 
that a critical stress, say, the yield stress, exists and that the load' stress is 
greater than the critical stress. Then after an almost instantaneous elastic 
yielding to the applied stress, the material will continue to deform in plastic 
deformation. Sketched schematically, in Figure 6-1, is the nonelastic or creep 
strain as a function of time. The creep strain is divided into three episodes, 
called the primary, the secondary, and the tertiary regimes. The tertiary part 
of the response is followed by brittle fracture (indicated by F in Figure 6-1). 

In the Earth this part of the strain history would be called an earthquake. Much 
of the present activity in the search for earthquake premonitors is connected 
with the local instrumentation of earthquake faults. A good example is the 
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San Andreas Fault where measurements of strain are being taken at this time to 
see if any accelerated creep is observable. If this is the case, then presumably, 
a region which is undergoing tertiary creep has been found; and brittle fracture, 
i. e. , an earthquake, is expected to be imminent on some undetermined time 
scale. 

But the problem, as it occurs frequently in nature, is hardly consistent 
with the simple explanation previously outlined, i. e. , with the premise that the 
observation of tertiary creep is the solution to the problem of earthquake pre- 
diction. The more detailed problem is rather simply illustrated by the fact 
that along some parts of the San Andreas Fault some creep is presently being 
observed while other parts of the fault are "locked," i. e. , at present no relative 
motion is currently being observed. 

Some questions present themselves. Why should one part of the San 
Andreas Fault be locked and not another? Is the locking due to the presence of 
stresses that are below the yield stress of the respective materials? Or is the 
material in the locked sections of the fault so durable that it does not deform 
significantly under even relatively large stresses? The problem then becomes 
one of determination of the shape of the response curve sketched in Figure 6-1. 

At the present time it is not known where a given region of the Earth is in the 
strain-stress -time space depicted in Figure 6-1. Further, the shape of the 
creep response curve itself, namely how long the various episodes of creep 
should last, is not known. 

The three curves of Figure 6-2 are intended to be schematic curves of 
creep response under constant load stress for a given material undergoing 
brittle fracture. For higher load stresses, the creep rates are higher, the 
time to fracture is shorter, and the strain at fracture is probably lower. 
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Figure 6-2. Strain Versus Time 


The situation is even more complicated when different materials are 
considered. For example, consider a given supercritical load stress applied 
to a material such as chewing gum or silly putty. The creep strains become 
large in short time, and fracture occurs relatively quickly when compared with 
geological situations (see Curve A of Figure 6-2). In the seismic case, for the 
same load stress discussed before, the response may depend seriously on 
frictional forces within the fault. For lower friction, as depicted by B in 
Figure 6-2, there may be some creep slippage, but as indicated before, it is 
not known whether this implies a short time to rupture as illustrated at point 
Bl, or a longer time as at B2. On the other hand, for a fault with high fric- 
tion as illustrated in Figure 6-2, the material may not deform significantly for 
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a long time until an episode of accelerated creep sets in, which may be very 
brief, as at Cl, or of a longer duration, as at C2. 

Curves such as those of Figure 6-2 also apply to parts of the Earth which 
are not on earthquake faults. The fact of the state of occurrence of an earth- 
fracture, i. e. , do not have earthquakes, implies that the critical stress is not 
reached. This may be due to the fact of occurrence of an earthquake nearby. 
Apparently knowledge of the state of strain in regions adjacent to earth- 
quake faults is important in understanding the state of strain in the faults 
themselves . 

Despite this rather bleak picture, some positive remarks can be made 
and tests can be proposed to ascertain the validity of various conjectures, and 
to clarify some confusing features of this problem. Multilate ration techniques 
cannot be expected to provide all the answers to earthquake problems, but cru- 
cial questions can be raised which multilate ration, by virtue of its unique capa- 
bility to obtain three dimensional station coordinates with precision, will no 
doubt be able to provide necessary inputs. 

6.2 PLATE TECTONICS 

According to the modern notions of plate tectonics, the surface of the 
Earth is covered by a relatively small number of rather rigid plates, on the 
order of 100 km thick. These plates* are presumed to be in relatively uniform 
motion with respect to one another, at least at places remote from their common 
boundaries. They are being created at the great mid-oceanic rises, by vir- 
tue of the efflux of mantle material at the rises; this material is cooled and 
attached to the edges of the plates as they recede from the rises. The plates 


*The Pacific plate, which is the largest of the plates, covers about 22% of the 
Earth's surface. 
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return matter to the mantle where the oceanic plates are in collision with 
continental plates or with other oceanic plates. At other plate boundaries, such 
as the San Andreas Fault which is the boundary between the Pacific and the 
North American Plates, the plates slide past one another. Here the mode of 
earthquake motion is consistent with strike-slip or horizontal motions of slip 
between two plates. At the trenches the earthquake motions have components 
consistent with compres sional stresses. 

Although the plates are presumed to be in rather constant relative motion 
in regions remote from the plate boundaries, this is not the case at the bound- 
aries themselves. Here the plates may be imagined to be held together by fric- 
tion, and motions take place abruptly at irregular intervals. These abrupt 
events are earthquakes and indeed a map of earthquake epicenters delineates 
the boundaries of the plates. According to the model of plate tectonics, the 
interiors of the plates are undoubtedly in a state of stress but these stresses 
are subcritical for fracture except in certain rare instances. Such an event, 
occurring in the interior of a plate, was the earthquake at Charleston, S. C. in 
1888. The earthquakes which occurred in Denver, Colorado, in 1962-1967 are 
testimony to the fact that the interior of the North American Plate is not stress 
free and indeed that the stress is subcritical at Denver; it was only when the 
critical stress for rupture was lowered by fluid injection at one point, that the 
series of earthquakes was initiated. It is further evident that the stresses in 
the interiors of plates are necessarily present in order that the plates can move 
relative to one another. 

Consider that the stress field in the neighborhood of a major fault can be 
represented by a model of a uniform block of the Earth's crust with a uniform 
strain-rate or stress field applied at infinity; the block has a fault within it and 
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sliding at the fault is prevented from occurring by frictional forces. 
Unfortunately the past history of earthquakes on the fault seriously creates an 
inhomogeneous stress, and hence strain field, in the vicinity of the fault. When 
a segment of the fault ruptures, stresses on the ruptured segment are lowered 
while stresses outside the torn segment are raised. The material in the neigh- 
borhood of a fault break probably undergoes some post-shock creep strain due 
to the change in the stress conditions; such post-shock creep has been observed. 
Hence astride a feature such as the San Andreas Fault, the stresses can hardly 
be considered to be uniform. As a crude simplification, there exists the prob- 
ability that following a great earthquake, there is a "waiting time" before another 
great shock will occur simply because the stresses and hence the strains in this 
region have been dropped. On this basis, a great earthquake should be expected 
to occur on that part of the fault which tore more recently in the past. 

From the above, it might be expected that a great earthquake might occur 
sooner along the San Andreas Fault of Southern California near its closest 
approach to Los Angeles instead of in the San Francisco area, since the last 
great earthquake in the former segment was in 1857 while in the latter case the 
most recent event was in 1906. However, these two regions are not in similar 
seismological settings. These dissimilarities in seismological settings have 
implications for differences in the frictional configurations as well, and hence 
for the respective creep histories. 

In reality, the contact between two plates such as the Pacific and North 
American plates in California is not a simple linear one -dimensional feature. 

The San Andreas Fault itself undergoes major changes in direction especially 
in Southern California. The Western United States and especially the seismic 
zones of California and Nevada are a rather complicated network of earthquake 
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faults, probably in some sense related with the plate boundary itself, or perhaps 

the network is the plate boundary. Many of the great earthquakes of California 
and Nevada history did not occur on the San Andreas Fault (Figure 6-3). A 
crucial question concerning the potential hazard on the locked (no observed 
creep) and the unlocked (proven or suspected creep) segments of faults concerns 
the sites at which the strains are being built up at this time. Estimates of the 
rates of motion of the Pacific and North American plates at places remote from 
their boundaries are of the order of 5 cm/year. In the neighborhood of the 
locked segments of the great faults, where is the strain rate of the order of 5 
cm/year if it is zero on the faults? In central California, on the "unlocked" 
segment of the San Andreas Fault, what is the significance of the fact that the 
slip rate is only about 1 cm/ year and where is the remainder of the strain field 
being absorbed? Does this imply strain build up on other important fault fea- 
tures of the Western U.S. ? The fundamental question to be answered is the 
extent to which the fields in the interiors of plates, and to the proximity of the 
plate boundaries, can be considered to be uniform. 

6. 3 QUESTIONS OF PRIMARY INTEREST 

A number of important questions have been raised in the discussion pre- 
sented in the previous sections. These include the following: 

1. Is the rate of relative motion of plate interiors truly uniform? Esti- 
mates of these rates are in the main derived from evidence which 
represent averages over several millions of years, although the 
recent JOIDES* deep-drilling results seem to imply that these aver- 
ages may hold for even longer periods; as long as 75 million years 


* Joint Oceanographic Institutions for Deep Earth Sampling 
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Figure 6-3. Large Earthquakes and Associated Faulting in the California 
and Nevada Region. The Locked and Unlocked Portions of 
the San Andreas Fault Zone are Delineated by 
the Shaded and Stippled Zones Respectively. 
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or more, But can these averages be applied to estimates of plate 
strain rates on a time scale of a few years? Measurements should 
be made of the relative motion of points on opposite sides of a plate 
boundary at points remote from the contact. It is suggested that mea- 
surement be made of the relative motion of the North American and 
European plates for this purpose. The Mid-Atlantic Ridge, which is 
the boundary in this case, appears to be less complicated than the 
Pacific-North American Boundary. Stations suggested for implemen- 
tation of a six-station experiment might include: New Orleans, 
Chicago and Halifax in the North American Plate and Rome, Green- 
wich and Stockholm on the European Plate. 

2. Are the motions of a given plate-interior relatively uniform, i. e. , 
does a plate-interior move as a relatively rigid body? Or are the 
stresses subcritical for fracture but still able to produce local creep 
in consonance with irregularities in the stress fields at the edges of 
the plates. In this case an experiment is proposed to see if there is 
any relative motion in the interior of the "stable" North American 
Plate. Station locations might include : Dallas, Kansas City, Minnea- 
polis, Atlanta, Ottawa and Halifax. 

3. If motions are relatively uniform in plate interiors, how close can a 
plate boundary be approached before these fields show a gradient? 
Here it is proposed to study the Pacific-North American Plate bound- 
ary, Proposed stations might include: Oahu, Farallon Islands, 

San Diego, Reno, Salt Lake City and Minneapolis. 
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4. What ground motions take place 100 to 200 km from a given fault 
zone? In this case local placement of six or more stations should 
satisfy all experiment objectives. Knowledge of the movements of 
distant points is a basic requirement in constructing a strain model of 
a major fault zone to be used as the basis of the understanding of 
earthquakes. Current geodetic-net techniques, while very accurate 
at short distances, are line- of- sight limited and are impractical 
beyond distances of about 100 km for the high accuracies needed. 
Multilateration, however, is not hindered in this respect. 

As an example consider how strain builds up along the San Andreas 
Fault zone as a function of whether the fault is locked or unlocked. 
Examination of Figure 6-3 shows four separate areas along the San 
Andreas Fault zone. These include: Northern California (site of the 
1906 earthquake), the unlocked segment to the south, the locked seg- 
ment farther south, and finally the Imperial Valley to the Gulf of 
California. Each of these segments will require different experiments 
due to the distinct geology peculiar to the respective areas. It is now 

believed that the type of slip that .relieve s_ strain . in the San Andreas 

fault zone varies along the length of the fault. In certain areas, some 
of the slip takes place as creep, that is, slow motion with little or no 
radiated seismic energy (earthquakes). In other areas, almost all 
of the slip occurs rapidly in large but infrequent earthquakes. The 
former areas are considered locked inasmuch as most of accumulated 
strain is being elastically stored for sudden release in further large 
earthquakes. The regional strain field along the San Andreas is 
strongly affected by whether the fault is locally locked or unlocked. 
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and any strain measurement program should take this into account. 

For example, in areas of fault creep, it is important to know if the 
input strain, measured by the displacement of points distant from the 
fault (100 to 200 km), is relieved by creep along or near the fault, 
which is relatively easy to measure. If the strain is relieved mainly 
by creep, then the likelihood of a large earthquake along that part of 
the fault is small. However, if only a small portion of the input strain 
is relieved by creep, then unlocked portions of the fault are also 
susceptible to large earthquakes. The San Andreas and associated 
fault zones in California and Mexico represent one of the most studied 
potentially destructive seismic zones in the world. Therefore, this 
region provides a truly unique location for getting relative station 
motions and, thus, regional strain data which are scientifically 
pertinent. 

5. What is the rate of strain build-up on other major fault zones? A 
number of catastrophic events have occurred in the Western United 
States other than on the San Andreas Fault. Note can be taken of the 
disastrous Owens Valley earthquake of 1872, which has been classi- 
fied as having a magnitude greater than the 1906 or 1857 earthquakes 
that occurred on the San Andreas fault. The rate of strain build-up 
on other major fault zones in the Western U. S. can and should be 
studied via the methodology of Multilateration, It is highly likely 
that the network of earthquake faults in the western U. S. will have to 
be considered interactively instead of individually. Here, the rela- 
tively low cost of Multilateration is a distinct advantage since instru- 
mentation of each of these faults for creep measurements may be 
extraordinarily complex. 


6-12 


JPL Technical Memorandum 33-605 



6. What is the driving mechanism for plate tectonic motions? It has 
been suggested that the failure of the Earth to have the equilibrium 
shape for a fluid rotating with the present angular velocity of the 
Earth is evidence for the presence of stresses within the mantle of 
the Earth. The estimates of the flattening of the Earth are 1/298. 3 
while the equilibrium flattening is 1/299. 7. This discrepancy of 
about one part in 200 is thought to be significant. 

The flattening is obtained from studies of perturbations in the orbits 
of satellites. A direct observation of the figure of the Earth may lead 
to better estimates of the stresses in the Earth. However, this 
information alone cannot indicate whether these stresses are caused 
by convection in the mantle or by finite strength of the mantle mate- 
rial. A direct observation of the figure of the Earth is possible by 
measuring the shape of the Earth at a geoid, or equipotential surface. 
Such an equipotential surface is provided by the surface of the ocean. 
Coastal stations at known heights above mean sea level could be used 
to determine the geoid. Suggested stations include: Halifax, Wash- 
ington, Puerto Rico, London, etc. 

6. 4 EARTH/ OCEAN TIDE CALIBRATION 

As mentioned in the previous sections, the prediction of earthquakes 
involves a succinct number of factors. One of these factors is the ability to 
make precision measurements across faults and between adjacent plates. 

These measurements are related to the plastic deformation of local ground 
areas. However, the Earth also experiences elastic deformation due to the 
perturbative effect of the Sun and Moon. Specifically, the acceleration of the 
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Sun and Moon on two points separated by a sizeable distance 1000-10000 km 
will have unequal components and thus will cause residual movements between 
the two points. Hence, as a function of time two points on the surface of the 
Earth will have a cyclic motion with amplitudes ranging from 1 to 30 cm. rela- 
tive to some mean datum. Furthermore, ocean tides will also cause differen- 
tial ground movements. 

It therefore follows that to obtain precise station coordinates the magni- 
tude of the Earth/ocean tides must be calibrated at each station. This calibra- 
tion, i. e. , the determination of curves of the sort displayed in Figure 6-4, does 
not present any serious problem with respect to obtaining precision plastic 
ground movements. The required calibration can be performed with the aid of 
a moveable gravimeter (placed at each station for a suitable period of time) or 
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Figure 6-4. Earth/Moon Tides 
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perhaps directly by analytic calibration. It is estimated that the elastic defor- 
mations in the crust of the Earth can be determined to the millimeter level. 

6. 5 CONCLUDING REMARKS 

This section has indicated the obvious: namely, that the prediction of 
earthquakes is a difficult task. It has been pointed out that Multilateration, 
along with carefully planned experiments will be a unique advantage to the pro- 
cess of earthquake prediction. Some preliminary station locations fundamental 
in answering the questions posed have been selected. Further analysis will 
permit the exact viability of these station locations to be established. 

It can be concluded, at least on a preliminary basis, that use of 
multilateration techniques will lead to a better understanding of earthquake 
mechanisms. Even if the consequences of the multilateration systems discussed 
herein lead to the ability of only roughly estimating the epoch of major earth- 
quakes, this will result in the development of long range planning to aid in pre- 
paring for such events. This long range planning, i. e. , setting in motion the 
political and financial plan for the predicted event will result in great financial 
savings because manpower levels need only be increased when the probability 
of an earthquake is high. Furthermore state officials could easily budget funds 
for use in future reconstruction. Short range planning, i. e. , the development 
of strategy for an earthquake predicted to occur within a week or perhaps days 
needs no explanation. 
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SECTION VII 


GENERAL THEORY OF MULTILATER ATION 
This Section is primarily concerned with the mathematical theory of 
multilateration, i. e. , the determination of three-dimensional station locations 
and orbital ephemeris by the use of geometric principles only. 

The first portion of this Section discusses the necessary conditions for 
the existence of a solution to the ranging problem. It is indicated that the 
sufficient conditions are not established but a preview of what will occur when 
systems comprised of 4, 5 and 6 stations are used is outlined with the aid of a 
schematic block diagram. 

In Paragraph 7.2 a derivation of the actual geometric equations is 
undertaken. It is shown that the solution to the ranging problem consists in the 
simultaneous determination of the intersection of a set of spheres each centered 
at the location of a station. These equations can be reduced at a fixed time to 
yield a single nonlinear equation which describes the overall geometric relation- 
ships. The repeated application of the reduced or fundamental equation at 
a repeated number of times subsequently will yield the solution of the multi- 
lateration problem. 

The third portion of this Section discusses the sensitivity of the 
solutions, i. e. , for fixed stations and orbital configurations, the manner in 
which errors in the observed range measurements map into the actual station 
positions is outlined. 

Paragraph 7. 4 introduces the reader to the degeneracies peculiar to 
the type of geometric schemes proposed herein. It is in this Section that the 
reasons for the degeneration of the 4 or 5 station configurations are explained 
in sufficient detail such that the interested reader can proceed to Appendix B 
for the complete theory, if it is so desired. The reason for the non-degeneracy 
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of the 6 station configuration also is stated. Finally, Paragraph 7. 5 explains 
why a solution to the multilateration problem is possible using three collinear 
stations. 

7. 1 NECESSARY CONDITIONS FOR THE SOLUTION OF 

RANGE-ONLY SYSTEMS 

As will be discussed presently, the mathematical relationships 
utilized in geometric ranging techniques give rise to one equation per station 
per strike of the satellite, which contains as parameters, the unknown position 
of the satellite [x,y, zj along with the unknown position of the station [X, Y, Z], 
and as known quantities, the ranges. Owing to the adopted coordinate system 
(see Section 4), in which the coordinates of the first three stations are given 
respectively by [0, 0, 0],[ X £ , 0, 0] , and [X 3 ,Y 3 ,0], the number of unknown 
station coordinates is reduced by 6. Hence, in the first three equations, only 
three unknown station coordinates are introduced. 


These facts are tabulated in Table 7-1. 

Table 7-1. 

Summary of Conditions 


Station Number 

Number of Satellite Coordinates 
Per Strike 

Number of 
Station Coordinates 

1 

3 

0 

2 

3 

1 

3 

3 

2 

4 

3 

3 

5 

3 

3 

6 

3 

3 
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In the previous table, the dashed line indicates that, after the third 
station, the number of station unknowns is always increased by three for each 
new station added to the system. These statements can be expressed mathe- 
matically as : 

UNKNOWN STATION COORDINATES = 3(1-3) +3 , I»3, (7.1.1) 

where I is the number of stations in the station net. 

The remaining unknowns, i. e. , the satellite coordinates, are simply 
given by three times the number of strikes (See Table 7-1) so that: 

UNKNOWN SATELLITE COORDINATES = 3N (7. 1. 2) 

where N is the total number of strikes. The total unknowns are therefore 
given by 


UNKNOWNS = 3(1-3) +3 +3N. (7.1.3) 

Now consider the number of measurements which are made at each 
strike of the system. At any given time each station makes one measurement 
of the slant range, p ; therefore the measurements define a sequence equal to 
(See Table 7-1): 

I, 21, 31 NI. 

Therefore, the total number of knowns is 

KNOWNS = NI , I>3, N = 1,2, 3, . . . (7.1.4) 

The knowns must be at least as many as the unknowns in order to yield a 
theoretically deterministic system, therefore: 

NI2 3(1-3) +3 +3N. (7. 1. 5) 
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Satisfaction of inequality (7. 1. 5) is a necessary condition for the 


mathematic solution of a system of equations determining the satellite coordi- 
nates at any given time together with the station coordinates. The solution of 
the inequality is given parametrically in Table 7-2 using I, the total number of 
stations, as the independent parameter. 


Table 7-2 

Parametric Solution of Inequality 


N = 31 ' 6 
" 1-5 


Minimum Total Number of 

Total Number of Stations (I) 

Strikes (N ) 


min 

3 

CO 

4 

6 

5 

5 

00 

vD 

4 

00 

3 


It should be emphasized that satisfaction of inequality (7. 1. 5) is only a necessary 
condition that the system under study be solvable^. It is only through detailed 
analysis, such as will be described in the following Sections and in Appendix 2 
that the sufficient cond itions peculiar to the ranging problem can be established. 

To provide the reader with a preview of the complex logic contained in 
geometric ranging schemes Figure 7-1 is introduced at this point. The 
figure outlines the fundamental degeneracies of the laser ranging problem and 
is divided into two branches depicting four /five, and six station ranging systems, 

*The previously introduced three station system is not strictly a range only 
system (See Paragraph 7. 5). 
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CONDITIONS WHICH 
CAUSE 

DEGENERACIES 



Figure 7-1. Degeneracies Causing System Malfunction 
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Block 1 indicates that there are certain orbital configurations in which 
the measured ranges do not change sufficiently and therefore whether the 
6 station (Block 2) or 4/5 station (Block 3) system is utilized - the mathematical 
equations will become degenerate and no solution will be possible. If the 
degeneracy does depend on the orbit (Block 4), then by using two circular orbits 
of different altitude or an elliptic orbit (Block 5) the degeneracy can be removed. 
The six station system, whether or not the conditions outlined in Blocks 6, 7 
and 8 exist will now yield a well defined numerical solution. 

On the other hand, the 4 or 5 station case requires that the question 
in Block 9 be answered. If the stations are nearly contained in a plane, 
then perhaps a new station sufficiently removed from the plane of the other 
three can be used in the station net. Similarly, if the distance between the 
stations is too short then intercontinental positioning of the stations will be 
required (Block 10). Satisfaction of these two conditions will also yield a 
workable system; otherwise the system using 4 or 5 stations becomes unwork- 
able regardless of the number of measurements which are taken. 

The previous rules of geometric ranging techniques are not obvious 
and the proof of the previous statements is involved. Instead of discussing 
these details at this point, it will be better for the reader to understand the 
basic fundamentals (Paragraph 7. 2, 7. 3) first and then, after reading the intro- 
duction in Paragraph 7. 4, proceed to Appendix B, if further details are desired. 
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DERIVATION AND METHOD OF SOLUTION OF MULTILATERATION 
EQUATIONS 


7. 2 


As will be seen presently, the derivation of the fundamental equations 
of the geometric method is very simple; however as explained in Appendix B, 
understanding the connection between the pertinent variables is a much more 
complicated matter. 

As illustrated in Figure 7-2 consider a station which makes a 
strike, i.e. , a measurement of the distance between itself and an orbiting 
satellite*, what geometry is described? Of course, the range distance 



Figure 7-2. A Station Measuring the Range Between Itself 
and an Orbiting Satellite 

*The word satellite is used throughout for simplicity; however, for certain con- 
figurations it should be remembered that an airplane would work equally as 
well as a satellite. 
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can be pictured as the radius of a sphere centered on the station and touching 
the satellite. Hence the fundamental equation peculiar to all geometric tech- 
niques, as depicted graphically in Figure 1-3, is that of a sphere which has 
been translated to a station located at coordinates X, Y, Z relative to the 
adopted origin at 0. Mathematically this relationship can be stated from the 
principles of analytic geometry as 

(x - X) 2 + (y - Y) 2 + (z - Z) 2 = p 2 , 
where x, y, z are the coordinates of the satellite. 

More specifically, let there be a total of I stations on the ground all 
of which range simultaneously to a total of N satellite trajectory points at 
different times, t. Then the following equations are self-evident, i.e. , 


(x 


n 


X.) 2 + (y - Y.) 2 + (z - Z.y 
1 *n l n l 




= 1 , 2 , ••• , 1 ) 
= 1,2. ••• ,Nj 


(7. 2.1) 


where 

til 

(X., Y., Z.) = the coordinates of the i — - station, 

l l l 

th 

(x^, y n> z^) = the coordinates of the n — trajectory point, 

th 

(P^ n ) = the slant range from the i— station to the 
n— trajectory point. 

It is convenient to work in the adopted geometric coordinate system introduced 

in Paragraph 4.4. In this coordinate system, the locations of stations 1, 2 and 3 

are used specifically to define the coordinate axes. Thus, by definition, 

Xj s Yj s Z 1 = Y 2 = Z 2 = Z 3 = 0. (7.2.2) 

The unknown station coordinates to be determined are X 1( X_, Y-, X ., Y., 

2 3 3 4’ 4 

Z 4 , - - - - Xj, Yj, Zj, a total of (31-6 ) parameters. 


7-8 


JPL Technical Memorandum 33-605 



For purposes of convenience and algebraic manipulation range square 
difference equations will be formed by subtracting Eq. (7. 2. 1) written for i = 1 
and for i = Z, 3, 4, * • ■ , I, namely, 


■ Zx X 0 + X" = 6, 
n Z Z Zn 


• Zx X, - Zy Y. + X? + yJ = 6, 
n 3 7 n 3 3 3 3n 


(7. Z. 3) 
(7. Z, 4) 


-Zx X. - Zy Y. - Zz Z. + X 2 + Y 2 + Z 2 = 6. , j = 4, 5, • • > , I , (7. Z. 5) 

n J n J n J J 3 J J n 

where the quantities 6^ n> (i = Z, 3, * * • ,1, n = 1, Z, • • • , N), are defined by 


, _ Z Z 

6. = p. - p. 

in in In 


(7. Z. 6 ) 


The differenced set of equations (7.Z. 3) through (7. Z. 5), together with 
the equation 


x 2 + y 2 + z 2 = P 2 n = 1, Z, 3, . . . N (7. 2. 7) 

n 7 n n In 

is completely equivalent to the original set of ranging equations (7, 2, 1), 

As mentioned previously, the advantage of using the range square 
difference equations is that the satellite variables can be eliminated conven- 
iently. Thus, from Eq. (7.2.3), x n can be obtained as 

X n = (X 2 - 6 2„»/ (2X 2>- < 7 - 2 - 8 ' 

Substitution of (7.2.8) into (7.2.4) allows y n also to be expressed as a function 
of the station coordinates and range data, i.e. , 


y 


n 







/(2Y 3 ). 


(7.2.9) 


Equations (7. 2. 7) can now be used to express as a function of the station 
coordinates X^, X^ and Y^: 
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(7.2. 10) 


z „ = [ P L - x n' X 2> - ^ (X 2’ X 3' V]'/' 

Note that z must be positive. A negative z n would correspond to a satellite 
trajectory point below the horizon from which a range measurement would not 
be possible. 

The expressions (7. 2. 8) through (7. 2. 10) can then be substituted into 
(7. 2. 5), giving a system of (I - 3) x N nonlinear equations in all the station 
coordinate variables: 

F (S) e X 2 + Y 2 + Z 2 - 2x (X,)X. - 2y (X_, X., YJY. 
m - j j j n 2 j *n 2 3 3 j 

- 2z (X,, X,, Y,)Z. -6. = 0, (7.2.11) 

n' 2’ 3 3 j jn 

where the index j runs from 4 to I, the index n from 1 to N, the index m from 
1 to (I - 3)N, and the vector S denotes the collection of all the unknown station 
coordinates (X^, X^, Y^ - - - - X^, Y^., Z^). 

The system of equations (7. 2. 11) can be solved easily by the Newton- 
Raphson method [7. 2]. More specifically, denoting a first estimate of any station 
parameter to be determined by S°, then an improved value of S^. is given by 
correcting S° by the amount AS^ obtained by solving the following system of 
linear equations: 


31-6 

s 


k= 1 


9F 

m 

3S k 


(S°)AS k 


= -F (S°), 

m — 


m = 1, 2, 


(I - 3)N. (7. 2. 11) 


The appropriate algorithm for this purpose is given in Appendix A. Numerical 
simulation has shown that the algorithm is extremely efficient, yielding accurate 
solutions in 1 or 2 iterations even from starting values which are 1 kilometer 
in error. The equations peculiar to these parametric schemes are simple, at 
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least from a conceptual point of view. The reader, however, is cautioned that 
the singularities inherent to these systems should be well understood prior to 
attempting to extract meaningful information. These singularities are discussed 
in Paragraph 7. 4 and in Appendix B. 

7. 3 SENSITIVITY ANALYSIS OF THE MULTISTATION SOLUTION 

In Paragraph 7. 2 it was shown that accurate solutions to the multilater- 
ation equations can be obtained in a straightforward manner. It now remains to 
analyze whether the solution is sensitive to small changes in the range measure- 
ments, This analysis is important because the measurements are subject to 
random and systematic errors which are unknown system variables. Speci- 
fically, if slight errors in these measurements produce large changes in the 
calculated station locations, then the mathematical method of solution and the 
overall system would not be useful in meeting the objectives of the present 
study, i.e. , in the determination of station locations accurate to the centimeter 
level. 

Numerical studies have shown that different configurations of the 
stations and the satellite trajectory points will in general result in different 
magnification factors (sensitivity numbers) in the error mapping from the 
measured data to the values of the station coordinates whose values are to be 
determined. In fact, there exist certain configurations which have infinite 
error factors, i.e. , geometric instabilities such that a unique solution of the 
multilate ration equations is not possible, even if a system of as many equations 
as unknown is formulated. Such cases are defined herein to be degenerate 
case s . Configurations of stations and trajectories that are close to being 
degenerate consequently will produce very large sensitivity numbe r s . Both 
degenerate and near-de gene rate configurations should be avoided in order for 
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the system to yield reasonable sensitivities, i.e., small error magnifications 
and thus small errors in the station coordinates. 

In order to decide whether or not a given configuration is geometrically 
stable, a sensitivity analysis must be performed. In principle, this type of 
analysis is simple although the algebraic details can become complex. Only 
the basic formulation of the analysis will be discussed herein. The details of 
the computational algorithm are discussed in Appendix A. 

Suppose each of the measured slant ranges, p^, is subject to an 
unknown range error, dp^. The resulting solution vector of the station coor- 
dinates, namely, 

S = (X 2 , X 3 , Y 3 , ••• ,X r Yj, Zj), (7.3.1) 

will then differ from the true solution by the amount 

dS = (dX 2 , dX 3 , dY 3 , ••• , dXj, dYj, dZj). (7.3.2) 

The purpose of performing a sensitivity analysis is to express dS> as a function 
of dp^ n and of the fixed geometry peculiar to a given geometric configuration. 
Since only small changes in range are to be mapped into relatively small 
changes in the station coordinates (for nondegenerate cases), analysis via 
differentials will be adequate. To start the error mapping process, Eq. (7. 2.11) 
is differentiated yielding 

(2X . - 2x )dX. + (2 Y. - 2y )dY. + (2Z. - 2z )dZ. - 2X.dx - 2Y.dy 
J n J J n J 3 j j n j y n 

- 2Z j dz n " d 6. n = 0, j = 4, 5, • •• ,1; n = 1,2, •• • ,N. (7. 3.3) 

The quantities dx^, dy^, and dz n can be obtained from Eqs. (7.2.8) through 
(7.2. 10) by direct differentiation. Then by substitution into Eq. (7. 3. 3), an 
equation of the following form is obtained: 
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A^dX 0 + B^dX, + C^dY, + a ( ^dX. + b< j) dY + c^dZ = d<r^\ 

n2n3n3njnjnj n 


j = 4,5, ••• ,I;n = 1,2, ••• , N, (7. 3. 4) 


where 


... / X.Y, - X.Y. 

' X A 1 X Z Y 3 \ 





(7. 3.6) 

(7. 3. 7) 

(7. 3. 8) 
(7. 3. 9) 
(7. 3. 10) 


/X, Y.X. x_ Z, y X,Z.' 


o-O) = I d6 ( — 1 _ ■) 3 _ — — 1 + — B . 3 _J 
n 2 2nyX z X.Y. z X. r z X.Y. 


X,Y, z X_ 
2 3 n 2 


z X-Y- 
n 2 3 


+ 2 d6 3n \ Y 


Y. 

A 

3 


y z. ' 

1 n 


Z. 


1 


2 

n 3 , 


z p ln dp ln " I d6 jn' (7 - 3 - U) 
n J 


For the purpose of simulation, dp^ n can be assumed to be a random 
variable with a Gaussian distribution having mean 0 and standard deviation of 
1 cm. The combined error vector dtr^ is evaluated according to Eq. (7. 3. 11) 
and the linear system, Eq. (7. 3.4), is solved to yield the corresponding station 
coordinate error vector dS. The components of dS divided by 1 cm are defined 
herein to be error-magnification coefficients. These coefficients will vary for 
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different sets of range errors even though the range errors are taken from the 
same distribution, i. e. , they are themselves random variables. The probabil- 
ity distributions for these coefficients, which are virtually impossible to derive 
analytically, except for normal random errors and stationary configurations, 
may be studied numerically by simulating a large number of identical configu- 
rations, i. e. , by repeated trails, with different sets of random range errors. 

A geometric configuration can then be shown to be stable only if both the mean, 
t, and the standard deviation, tr, of each error magnification coefficient are 
small numbers consistent with the accuracy requirements of the present study 1 *'. 
Typical study results have been presented in Section V. 

It should be noted that the above procedure is possible only if the par- 
ticular system configuration under study is nondegenerate. For degenerate 
situations, the coefficient matrix corresponding to the linear system, i. e. , 

Eqs. (7. 3.4), will be rank- deficient. This means it will not be possible to 
solve for the station error vector dS uniquely. A number of degeneracies have 
been discovered in the present study. These cases (see Paragraph 7. 1) are dis- 
cussed further in Paragraph 7.4. However, the most important degeneracy 
deserves to be mentioned in this Section. This is the coplanar degeneracy 
associated with the 4 station configuration. 

As discussed in Paragraph 7. 1, four stations with a minimum of six 
range measurements provide, at least in principle , enough equations for all the 
unknowns (station as well as satellite coordinates) to be determined. However, 

^Specifically, if the indicated calculations are performed M times using differ- 
ent sets of random range errors (drawn from the same distribution), with the 
mth calculation producing a set of sensitivity coefficients ds( m ) then 
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if all the stations are in a plane, then the system is degenerate no matter how 
many measurements are utilized for the solution. In fact, this degeneracy has 
forced the originally proposed 4- station case to be abandoned as a system for 
short baselines because stations are nearly in a plane unless they are separated 
by intercontinental distances. Only if at least 6 stations are used can this 
degeneracy be overcome. It should be mentioned that the 4-station system is 
still adequate if used intercontinentally. The next Section is intended to pre- 
sent a more detailed discussion of the problem. 

7.4 FUNDAMENTAL CAUSE AND REMOVAL OF SYSTEM 

DEGENERACIES 

In this section the fundamental causes for some of the degeneracies 
previewed in Paragraph 7. 1 are discussed in greater detail; i.e. , the main 
reason for the degeneracy due to the number of stations is isolated. The 
detailed discussion of this problem will be found in Appendix B. The reader 
should consider this section as an outline which will serve to aid in the further 
understanding of this problem. 

It will be demonstrated herein that the 4 and 5 station systems are 
always degenerate for geometries wherein the stations are contained in a plane, 
whereas the 6 station system is generally well defined. 

Consider the resolvent equation of the multilate ration process, i. e. , 
Eq. (7.2. 11) written for a system of 4 stations. For purposes of clarity the 
equation will be written out explicitly for one given time point or strike. The 
result is displayed in Figure 7-3. 

Examination of Figure 7-3 immediately yields the Second Solution 
(Section I) of geometric ranging, that is, repeating the resolvent equation at 6 
or more different time points will result in an algebraic system of six or more 
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Figure 7-3. Fundamental Equation 





equations in the six unknown station coordinates. Theoretically, it is evident 
that a solution to this problem exists. 

A cursory examination of the fundamental equation displayed in Fig- 
ure 7.4. 1 shows that if = 0, the equation collapses to the following form: 

a + 6 2 (3 + 63 Y = & 4 . (7. 4. 1) 

where a, (3, and Y are strict functions of the station coordinates. (The explicit 
formulas for these quantities are given in Appendix B. ) Equation (7. 4. 1) is 
obviously a linear equation with known coefficients 63 , and unknown param- 

eters a, (3, and Y. Let three strikes of the satellite be taken, i. e. , let Equa- 
tion (7. 4. 1) be written at three separate time points. The result is 

* + & 2 n P+ 6 3n Y = 6 4n’ n = 1, 2, 3, (7.4.2) 

a system of three equations in three unknowns which always has a unique 
solution unless the determinant of the system vanishes. The determinant of 
the system, i. e. , 



1 

6 21 

*31 1 

| 


D = 

1 

6 22 

6 32 

(7.4. 3) 


1 

6 23 

6 3 3 



is a pure function of the range data and could conceivably vanish for certain 
combinations of the data but this occurrence would be coincidental, i. e. , due 
to peculiar combinations of station and satellite coordinates and can be ignored 

*A discussion of this kind of degeneracy will be found in Appendix B. 
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herein. Hence under the assumption that D 4 0, system (7. 4. 2) can always 
be solved but only three quantities can be obtained, namely, the parameters 
a, (3 and Y . Taking a new strike on the satellite only yields a new equation of 
the same form as (7.4. 2) and results in an overdetermined system of linear 
equations. Obviously a least squares solution can then be obtained to the sys- 
tem but only for the three unknowns. Hence, no matter how many data points 
(strikes) are processed only three quantities can be obtained. However, when 
Z^ = 0, there are five unknown station coordinates and therefore a unique solu- 
tion for the explicit coordinates is not possible. 

Consider the addition of a fifth station; obviously, the equation displayed 
in Figure 7. 4. 1 can be written with X^, Y^, replaced by X,., Yg, Z^ and 
in the planar case if Z^ = 0 another independent linear system such as 
Eq. (7. 4. 2) can be obtained. This independent linear system upon solution 
will determine another three unknowns distinct from the parameters previously 
determined using the first four stations . But in a system of five stations 
(planar) there are seven station coordinates, and as discussed above only six 
independent unknowns can be obtained. Obviously a solution is not possible. 

However, when the sixth station is introduced into the system, then by 
the previous reasoning another three unknowns distinct from the first two sets 
of unknowns will be determined, and since in the planar geometry only the two 
additional unknowns [X^, Y^] are introduced, it follows that: nine independent 
coefficients of the coordinates have been specified for the nine unknown coordi- 
nates so that unless coincidental degeneracies of the kind discussed in 
Appendix B are present, a unique solution for the station coordinates now exists . 

If the Z^ coordinates are not exactly zero, then the equation displayed 
in Figure 7-3, written at six different times will permit a unique solution for 


all the coordinates; however unless the Z^ component is large enough a nearly 
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singular numerical system will be obtained. This should be evident since if 
= 0 the preceding results show that the analytic solution degenerates 
absolutely. In the five station case the same result can be inferred since the 
mathematical degeneracy is still present in the equations. However when the 
sixth station is added the degeneracy is removed generally and then, even for 
small Z^, Zg, Z^, a well conditioned solution exists. 

The coplanar degeneracy and other forms of degeneracy which can be 
readily overcome are discussed further in Appendix B. 

7.5 THREE-STATION SOLUTION: A DEGENERATE BUT SOLVABLE 

CONFIGURATION 

In Paragraph 7. 1 it was proven that at least four stations are neces- 
sary for a mathematically deterministic solution to a truly range-only system. 
However, if additional information is available about station configurations, 
then less than four stations conceivably may be adequate to form a workable 
system. Indeed, it will be shown presently that the fundamental multilate ration 
equation (see Fig. 7-3), written for three of the four stations aligned in a 
straight line, degenerates into a form such that the baselines between the col- 
linear stations can be uniquely determined. 

To obtain the three-station solution from the general multilateration 
equation displayed in Fig. 7-3, consider the following substitutions: 

Y 4 = Z 4 = 0 (7. 5. 1) 

Physically, this means that Station 4 is placed on the X-axis of the adopted 
coordinate system, i. e. , Stations 1, 2, and 4 are now collinear. The result 
of this substitution at any fixed time is 

X 

X 4 -( X 2‘ 5 2 ) 4' S 4 = °- < 7 - 5 - 2 ' 
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which can be written as 


X 4 6 2 ‘ X 2 6 4 + X 2 X 4 (X 4 " X 2 } = °‘ (7> 5 ‘ 3) 

Equation (7. 5. 3), written at two separate times, will provide two 
independent equations for the two unknown baselines X-, and X^. Hence, a 
unique solution can be obtained by taking two strikes. Note that the coordinates, 
X^ and Y^, as well as the pseudo measurement 6^ associated with Station 3, do 
not enter into Eq. (7. 5. 3). This in effect means that only three stations (1, 2, 
and 4) need be employed for the solution. Furthermore, it should be noted that 
only baselines are determined by this solution. 

An alternate approach to the three - station solution independent of gen- 
eral multilateration theory will be discussed in Appendix C, where more mathe- 
matical details are given. 
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SECTION VIII 


OPERATIONAL CONSIDERATIONS 

ffinm)) This Section primarily deals with the topics which impose constraints 
on system operation when the basic principles of the multilater ation technique 
are implemented in real time. 

Some of the highlights of this Section include: the manner in which 
orbital configurations must be selected in order to achieve a numerically 
stable solution to the multilater ation equations; the method of smoothing the 
data by use of a reference orbit; an examination of the accuracy to which 
station clocks must be synchronized; and the discussion of how the trajectory- 
information required for aiming of the laser unit can be extracted as a 
by-product of the solution of the multilateration equations. The methods for 
selecting the mutual viewing windows of the system stations, correcting for 
the effects of atmospheric attenuation on the laser beam, and the examination 
of laser safety considerations relative to humans, are analyzed, and pertinent 
recommendations are given. 
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8. 1 ORBIT ELEMENT SELECTION 


Any experiment utilizing a satellite whose fundamental objective will be 
to map the relative station locations of a set of ground stations via three- 
dimensional laser fixes, will be fundamentally linked to the dynamic and 
orientation elements of the adopted satellite orbit. 

8. 1. 1 Dynamic Element Selection 

The shape, or dynamic orbital parameters, such as eccentricity and 
semi-major axis, are important parameters to select because of the geometry 
constraints peculiar to any scheme which intends to invoke geometric principles 
to aid in the process of station location estimation. 

Numerical studies have indicated that for purposes of multilateration 
the selection of the orbit containing the satellite/r etr or eflector is of paramount 
importance. Even though the mathematical process of solution does not depend 
explicitly on the assumed orbit, it has been observed that the actual measure- 
ments taken over a viewing arc, must have a suitable dispersion in their 
respective magnitudes, i. e. , the slant ranges must undergo substantial change 
within the common viewing window (see Paragraph 8. 2) available to the station 
net. In fact, it has been observed that if a circular orbit is assumed, the 
multilateration equations yield ill defined estimates for the station coordinates. 
However, if a system of two circular orbits with well separated radii is used 
this problem is overcome and subsequently the system of equations becomes 
strongly non- singular . Hence, one solution to the orbit selection process is 
to adopt a system of two circular orbits, e. g. , with altitudes of 500 and 75 0 km, 
and range over discrete segments of the viewing arc to these satellites. It 
should be noted that the satellites need not be visible to the stations simulta- 
neously. Operationally this implies that a batch of data would be processed 
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from one satellite first and subsequently the data taken by ranging to the second 
satellite would be processed and combined with the first batch of data. 

A second solution to avoid ill conditioning of the equations due to 
selection of a single circular orbit is to use an elliptic orbit having effective 
altitudes ranging from 500 to 750 km. This type of orbit, due to the combined 
rotation of the Earth and orbital perturbations acting on the orientation plane 
of the orbit (See Paragraph 8. 1.2), will have the same stabilizing effect upon the 
system equations as the use of two circular orbits. 

For purposes of yielding well defined geometries, it is also beneficial 
to select orbital altitudes which are commensurate with the baseline separations 
between stations. 

In summary, an operational system requirement will be the use of: 

• two circular orbits, or 

• an elliptic orbit, both with 

• altitudes commensurate with the relative station baseline separations. 

These operational constraints can be attained with relative ease. 

8. 1. 2 Orientation Element Selection 

The orientation elements, namely, the inclination of the orbit to the 
equator, i, the longitude of the ascending node, Q , and the argument of 
perigee, w, will now be briefly discussed (see Figure 8-1). 
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Figure 8-1. Orientation Elements 
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Since, the second zonal harmonic, J^, causes the elements and u> to 
have time rates of change defined by the following relationships' [4. 5] 


u> +/3 
0 2 


2 ,, 2 . 

a (1 - e ) 


, 5 . 2 . 

c - ^ sin l 


n (t - t Q ) 


n = n. 


2 /i 2 i 

a (1 - e ) 


~2 cos i ) n(t - t Q ), 


where and are the injection values of u) and , respectively, and the 
anomalistic mean motion n~ is defined as 


n 


k 


\Zp a 


-3/2 



3 

2 


. 2 
sin 



with 

3/2 

k = Gaussian planetary constant = . 07436574 (earth radii) /min, 
p = Sum of masses of satellite and Earth in terms of most ponderous mass, 
tp = Adopted epoch time, 

a systematic study of the desired injection values of ^2 q must be performed. 
Therefore by selecting the proper value of for a specified value of i, it will 
be possible to make certain that: first, the orbital track will pass directly 
over the optimum satellite/station geometric configuration, and second, that 
a convenient time duration between satellite injection and experiment initiation 
will be available to meet ground operation constraints. 


The effect of j£ on the orbital parameters, although small, is not negligible 
in orbital operations. 
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It must be remembered that the satellite track will only pass in exactly 
the optimum position over the stations one time per satellite/station repeat 
period. These synodic repeat periods are of the order of years. Assuming 
that was chosen such that the optimum geometric configuration would occur 
on revolution N, then, the geometric accuracy picture for multilater ation 
improves as the satellite approaches revolution N, is a maximum at revolu- 
tion N and slowly degrades to zero (no visibility) on a revolution by revolution 
basis after revolution N. When the rise/set geometry again becomes common 
to all stations, the same cyclic behavior of the accuracy will occur, however 
due to the resulting phase diffference between the satellite and Earth the 
accuracy of the geometry will be decreased until, as mentioned previously, 
the ground track precisely passes over the geometrically optimum station 
configuration, i. e. , until satisfaction of the synodic period is achieved. 

The final parameter, the inclination, must be selected such that 

i > <i>> 

specifically, such that the orbital inclination exceeds the most northerly or 
southerly station latitude, $ , by a sufficient amount to ensure geometric stability. 
Intuitively it would seem that a near polar orbital system would offer the advan- 
tage of covering all available ground areas whose location might want to be 
surveyed. Furthermore, the polar orbit causes the rate of change of the 
longitude of the ascending node to vanish (ft = 0) and thus simplifies the 
orbital injection proces s . 
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8.2 


VIEWING WINDOWS 


From an operational point of view, the supporting software to any 
multilate ration experiment must have the ability to define the appropriate 
viewing window within which the measured data, in this case time pulse or 
range data, will be common to the system stations and the satellite at any given 
time. 

Figure 8-2 illustrates, in cartoon fashion, the rise/set histories 
of a satellite which can be expected from a net of four stations; obviously, a 
similar diagram can be displayed for a six station system. 



0 90 180 270 360 

ECCENTRIC ANOMALY (E), deg 

Figure 8-2. Superposition of Rise/Set Functions 
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The Figure is a graph of the rise/ set function, F, developed in [4. 5] 
as a function of satellite eccentric anomaly, E, i. e. , 


F = a(cos E - e) P • Z + a vl- e^ sin E Q • Z - G = 0, (8.2.1) 


where 


Z e cos cb c 
x 


Z = cos cb sin (0^ + 

y \ 0 


( e o + if ♦ t - t 0 j; 

d0 [E - e sin E „ 
dt L n 0 ] 


Z e sin cb 
z 


G = G^ cos^ 4> + G^ sin^ cj) 


with G j and G^ defined in Appendix J, and where for convenience the notation 
is repeated: 


4 > = 

H = 

6 0 = 
h = 

T = 

d0/dt = 

f = 


a 

e 


geodetic latitude of station, 

station elevation above and normal to adopted ellipsoid, 
epoch sidereal time at injection, 
anomalistic mean motion, 

time of perifocal passage or nodal crossing, 
sidereal rate of change, 
geometrical flattening of Earth 
equatorial radius of Earth, 


Note that the eccentric anomaly and time are simply related by Kepler's 
equation. 
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P = unit vector pointing from dynamical center to orbital perigee 
or nodal crossing, 

Q = unit vector in orbital plane advanced to P by a right angle 
in direction of motion. 

As can be seen from the figure the rise/ set function is a smoothly 
varying analytic function that achieves different maximums and minimums at 
different times for each station of the station net. The first zero of F^ repre- 
sents the rise of the satellite with respect to station i. From Figure 8-2 
it is simple to extract the rule for the limiting view period wherein all stations 
can see the satellite, namely, 

At = t* - t* > (8.2.2) 

where, if I is the total number of stations. 


* 

's = 

MIN (t sr t s 2 ’ • • • 

• ■ ’ w 

t* = 

R 

max ( t ^ , t 2 f • • 

• • • ’ ^1*' 


that is, t^ is the critical time at which the ranging experiment must start, and 
tg denotes the time at which the ground equipment can cease operation on a per 
revolution basis. 

The difference At is the duration of common visibility or system 
up- time on a given orbital pass. Obviously, due to the rotation of the Earth 
and the previously discussed orbital perturbations, the value of At will vary 
from zero to a maximum and then back to zero. The maximum value of At will 
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correspond approximately to the satellite overpassing the centroid of the station 
net. The time t^ corresponding to will define the optimum viewing 

window available for laser ranging experiments. 

Future predictions of these times present no difficulty and can be 
handled directly in the operational software package. 
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DATA SMOOTHING 


8. 3 

During actual operation of the laser ranging system a continuous 
stream of simultaneous fixes will be obtained from the set of ground stations. 

The laser measurement rate will be of the order of one measurement every 
six seconds. As is common in any measurement technique a given pulse will 
be corrupted with random noise and therefore occasionally bad measurements 
will be obtained. Since the hardware or measurement equipment will be cali- 
brated as to expected accuracy, i. e. , a standard deviation of the time pulse 
or range data peculiar to the equipment will be available, it will be a simple 
task to reject all obviously bad points by culling out any measurements whose 
residual is a factor of three greater than the standard deviation associated 
with the equipment accuracy' 1 '. Once this gross filtering process has been 
accomplished smoothing of the data can be undertaken. 

Two basic requirements must be satisfied when smoothing is to be 
performed. First, an analytically smooth function which represents the physi- 
cal behavior of the data must be adopted. Second, an analytic function which 
represents the data over a suitable duration with respect to the independent vari- 
able should also be assumed. 

The independent variable in three dimensional laser ranging will be the 
time, t, in minutes after some conveniently adopted epoch, t Q , say, midnight of 
the day when the experiment is to be performed. The accuracy to which t must 
be known, i. e. , to which the epoch clock synchronization must be specified is 
discussed in Paragraph 8.4. 


This is known in statistics as the "3cr criterion, " and for a normal dispursion 
corresponds to a confidence level of 99. 9 percent. 
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The first analytic function for data representation (it will be assumed 
that the time pulse data has been converted to actual range, p, or extension 
parameters) which comes to mind is of the form: 


Pi = a i + b.(t - t^) + c.(t - t^) + . . . , 


(8. 3. 1 ; 


where as discussed in Paragraph 8. 2, t is the maximum rise time associated 

-K 

with the four observing stations, and a., b., c^ are suitable least squares coef- 
ficients obtained, probably, with equal weights for all stations. This type of 
representation of the data is nothing more than a Taylor expansion about t^ and 
for an assumed number of terms will obviously degrade as a function of time. 
The previous function is excellent for synchronous satellites but is of limited 
use for low Earth satellites. 

To circumvent the previous problem and obtain greater accuracy in 
the number of significant digits consider the following scheme. 

Let a two-body or Keplerian orbit including secular perturbations (see 
Section 8. 1) be adopted as a baseline reference. To determine the elements of 
the baseline orbit to fair approximation, a very short arc of data can be fit as 
previously described via a Taylor expansion at three stations such that 


p. = a. + b.(t - t^) + c.(t - t^) 2 , i = 1,2,3, (8. 3. 2) 

and by differentiation 

p. = b. + 2c. (t - tT ) . (8. 3. 3) 

l l i R 

For purposes of discussion only quadratic terms are retained. 

Using a trilateration technique [l. l] , [4. 7] , it now follows that the 

six data measurements, i. e. , p^, jx at the middle of the short arc can be 
mapped into the elements a, e, i, Q , w , T , and thus the elements of the 
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baseline orbit become determined. These elements are now adopted as a 
base or reference . Being in possession of these elements permits the orbital 
evolution on a secularly perturbed Keplerian orbit [4. 5] to be performed and the 
range data to be computed from the reference orbit over the duration wherein the 
true range data is to be measured by each observing station. The accuracy of 
the Keplerian reference orbit is not of any importance for the scheme to be dis- 
cussed because the reference orbit ephemeris will be added in and then removed 
without corrupting the actual data. 

For all measurement times, n, at each station it follows that the fol- 
lowing residuals can be computed 


A 


P. 

m 



i = 1 , 2, . . . . I, 
n= 1,2 N, 


(8. 3.4) 


i. e. , the difference between the true measurement p and the computed Keplerian 
approximation p is formed on a point basis (sequentially). 

The recommended procedure is to fit the residuals A p as a function 
of time via least squares as 


Ap.(t) 


a. 

i 


+ Ta . (t 

l 


tR> + 


C i (t - 


t"' ) 

Tl ; 


i = 1 , 2, ... I, 


(8. 3. 5) 


where all calculations are carried throughout to 18 digit accuracy. Once the 
residuals are known as a function of time, the actual model for the range history 
as a function of time from each station can therefore be represented as 


P.(t) = P^(t) + Ap.(t) = A.(t) + b.(t - t*) + c £ (t - t*)^ 


(8. 3. 6) 


The use of a second order polynomial in the fit has been found to be adequate. 
See also the discussion at the end of the Section. 
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where 


A. (t) = (t) + a. • (8. 3. 7) 

It should be noted that in Equation (8. 3. 6), p^(t) is calculated by adding the 
previously derived Keplerian range to the fitted residual and thus the actual 
measurement is not corrupted by the reference orbit. 

An even stronger representation can be obtained by letting the curvature 
of all arcs be modeled simultaneously by least squares fitting of the data to the 
more generalized form 

Ap. = A + B (t - £)“ + C.(t - t*)P (8.3.8) 

1 1 1 K IK. 

where a and pare exponents determined via the fitting process (see Appendix G). 
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SIMULTANEOUS DATA DETERMINATION 


8. 4 

It has been shown in Section VII that the relative position of stations 
located on the surface of the Earth can be determined accurately if a sequence 
of slant ranges from a satellite is known accurately and simultaneously. For 
purposes of multilateration, it is important that the ranges from all stations 
to the satellite be known at the same time . But the distances between the 
satellite and any two stations differ at any given time, and since the satellite 
is moving during the transit time of a ranging pulse, it might seem difficult to 
achieve the necessary simultaneity of the range measurements. To overcome 
this difficulty one must use an interpolation scheme. Roughly speaking, if the 
range history as a function of time can be modeled, then slant range histories 
such as those displayed in Figure 8-3 can be generated. 



Figure 8-3. Range History as a Function of Time 
Therefore, by adopting some time t and entering the curves as in 
Fig. 8-3, a set of pseudo observations can be obtained which must be synchro- 
nized. In this section the relative accuracy of the curves which are required 
to obtain the pseudo observations will be discussed, i. e. , the accuracy of the 
epoch time peculiar to each station and therefore to the data will be specified. 


It will be shown that these requirements are not stringent. 
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Many transit times T(L) of a pulse emitted at times L are measured 
from each station independently during a pass, i. e. , during the viewing period 
common to all stations. Typically, for a pass lasting 15 minutes, a pulse can 
be fired every 6 seconds with existing equipment, giving 150 values for T(F). 

Let the range p be known to within 10 and 100 m from an a priori know- 
ledge of station location and simple trilateration. From the nominal range 
p^(t), it is possible to obtain the residuals : 

A(t.) = | T(t.) - p N (t ). (8. 4. 1) 

Now, with the data available it is possible to perform a least squares 
analysis to smooth the measurements (8. 4. 1). Accordingly, having selected a 
set of polynomials, i. e. , Chebychev, powers of t, or the like, for modeling the 
data, it is possible to write 

A(t.) = 2a n p n (t) = P(t), (8.4.2) 

and to determine the coefficients a^ by the method of least squares. In this 
manner a smooth time function P(t) is found which represents the residuals 
A(L). Naturally, P(t) is determined in this manner for each station separately. 
The determination of the range p(t) as a function of time from P(t) for each 
station then proceeds as follows. First, the relationship between the range 
p(t) and the total observed transit time T(t) must be known. It can be shown 
(Appendix E) that within the accuracy of a centimeter this connection is given 
by: 

T(t) = | p(t)[l + c' 1 p(t)] . (8.4.3) 


The use of the nominal range is not essential but is convenient for data 
processing. 
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In Eq. (8. 4. 3), c is the velocity of light and p the range rate. It must 
be emphasized that t in Eq, (8. 4. 3) is the time of emission of the pulse possess- 
ing the transit time T. 

It is now an easy matter to extract the range information from P(t). 

In fact, having substituted (8. 4. 3) and (8. 4. 1) into (8. 4. 2), differentiation of 
(8. 4. 2) reveals that 

P(t) = p(t) + c' 1 p(t ) 2 + c' 1 p(t)p(t) - P N (t) (8. 4. 4) 

so that 

P(t) = P(t) + p N (t) - j (p 2 + pp) . (8. 4. 5) 

Expression (8. 4. 5) is inserted into Eq. (8. 4. 3) to obtain the actual range from 
the measurements. Further analysis (Appendix E) will show that the last term 
on the right hand side of Eq. (8. 4. 5) is totally negligible. From Eqs. (8. 4. 3), 

(8. 4. 1) and (8. 4. 5) the range at time t is then given by: 

p(t) = [P(t) + p N (t)J [l + c _1 [P(t) + P N (t)]J . (8. 4. 6) 

The above expression contains only known quantities, namely, the nominal 
orbit range, and its time derivative, and the least squares fit P(t) of the data 
and its time derivative. If the transit time T is measured to an accuracy of 
0. 06 nsec, Eq. (8. 4. 6) gives the range between a station and the satellite 
within an accuracy of 1 cm. 

The procedure outlined so far is applicable to any one station. What 
is needed is a knowledge of the ranges from different stations simultaneously . 
This is tantamount to demanding that the station clocks be synchronized, i. e. , 
if the epoch times measured at different stations are the same, simultaneity 
is assured. 
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Suppose that two unsynchronized stations (1 and 2) measure the range 
to the same satellite. Let the clock at station 2 read the time t^ = t^ + 6 + a t^, 
where t^ is the time at station 1 and where 6 is the time difference (clock offset), 
and the difference in rate of clock 2 from clock 1, The parameters 6 and a 
are usually unknown, so that a comparison between and which is 

needed for achieving simultaneity, is difficult. However, as an approximation 

P 2 ( t 2 ) = P 2 ( t 1 ) + + at i )p( t 1 )’ (8.4.7) 

for small time differences t^ - t^. In Appendix E it is shown that p is at most 
about 1. 5 km sec ^ . Therefore, if 

(6 + Q?t) x 1. 5 (km sec ^) < 0. 5 cm (8. 4. 8) 


or if 


6 + at < 3 psec, (8. 4. 9) 

sufficient synchronization is assured. Quartz oscillators in conjunction with 
a ground wave LORAN system (See Glossary) are sufficient to achieve the 
desired synchronization, particularly since synchronization is only needed 
during an individual pass of comparatively short duration. Subsequent passes 
may be synchronized independently. 

In cases where a ground wave LORAN system is unavailable because of 
the long distance between a laser station and the nearest LORAN station, a 
satellite equipped with a radio beacon may be used. An S-band pulse trans- 
mitted to, for instance, two different stations will be received with a time lag 
given by: 


At = 



(8. 4. 10) 


In the above equation At is an exact measurement but the ranges are known only 
approximately due to errors in station location, errors in satellite orbit, and 
errors in the transmission medium. In order to stay within the allowed 3 psec 
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synchronization error, the errors previously listed must total no more than 
0. 9 km in equivalent range. The uncalibrated transmission medium errors 
(troposphere and ionosphere) amount to less than 20 m for elevation angles 
greater than 10 degrees. The errors in station location amount to no more 
than 30 m. This leaves an allowable 850 m for the indeterminateness of the 
satellite position, and it is not diffucult to have knowledge of the satellite's 
position to that accuracy. 

Therefore, it has been shown that a data smoothing procedure together 
with synchronization of station clocks within several psec ensures the necessary 
range simultaneity. 
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8.5 ATMOSPHERIC LASER CORRECTIONS 


This section discusses the method of correcting laser measurements due 

to atmospheric attenuation. These corrections are important when station 

positions are to be determined to centimeter accuracy. 

Laser beams are different from radio waves and are refracted in the 

atmosphere strictly due to dry air effects. This important fact is due to 

5 

relative frequency levels, i.e., the frequency of laser light is 10 times that 

of microwaves used for radio tracking and ranging. In this case the dipole 

moments of the water molecules respond negligibly to the radiation and hence 

14 

may be ignored. There are no other spectral lines in the neighborhood of 10 Hz 

(ruby laser light). As will be shown herein it will be possible to obtain accurate 

range calibration [8. 1] , [8. 2] in actual system operation. 

The dry zenith range correction, Ap , can be computed from either 

z 

surface measurement or radiosonde balloon data (See Glossary). For a static 
equilibrium atmosphere, the zenith range correction can be derived analytically 
as a function of surface pressure only, as [8. l] : 



11 . 6 R 
g 



(8. 5. 1) 


where Ap is in meters, R is the gas constant, P_ is the surface pressure in 

Zt u 

millibars, and in these units: 


g_ 

R 


34. 1 


K 

km 


Equation (8. 5. 1) is obtained under the assumptions of static equilibrium, perfect 
gas, and constant gravitation acceleration, g. If these assumptions are 
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adopted, the zenith range can be predicted with an uncertainty of ±2 mm for a 
precision of ±1 millibar in pressure measurement. Hopfield (after comparisons 
with observed data) [8. l] claims 2 mm as the uncertainty of the zenith range 
measurement from Eq. (8. 5. 1). It seems however, that the 2 mm uncertainty 
is too optimistic. First, the air convection or vertical current which is common 
especially near a lake or in the hot desert, may appreciably violate the assump- 
tion of static equilibrium. Second, there is no observed zenith range data to 
the above accuracy with which to perform a comparison to some standard. 
Observed data used by Hopfield must be computed from radiosonde balloon 
measurements which have a 3 cm uncertainty in zenith range correction. 

It has been shown [8. 6] that tropospheric range corrections can be 
expressed conveniently by the following equation: 

A P = cyj^ F (r , 0 o (r,Y))|l - ^ dr, (8. 5.2) 


where: 


R 

y 


1 + <»F 


e 

o 


(r,Y) 


radius of Earth, 

elevation angle of the ray path between the laser station and the 
satellite , 

refractive index of the troposphere, 

- Y + cos ^ - C ° S — -j = raypath expressed in spherical 

coordinates. 


_4 

ot = 3 • 10 (typically). 

F = scale factor depending on the modeling of the troposphere. 
Expression (8. 5. 2) has been used to obtain simplified but accurate formulas for 
the atmospheric range corrections. 
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The range correction, Ap, for a spherically symmetric atmosphere in static 
equilibrium, as suming no bending of the raypath, is relatedto the zenith value by 8. 3 

A P. 


Ap =■ 


sin y 


(8. 5. 3) 


where is the elevation angle, i. e. , the angle between the horizon and the line 
of sight to the satellite. 

However, bending effects should not be neglected for deviations from 
angles lower than 45 degrees. Equation (8. 5. 3) has an error of 1. 5 cm at 
= 30°, and 5. 3 cm at = 2.0°. For more accurate mapping valid to lower 
elevation angles, the following modified formula is recommended [8. 3]: 

A P„ 


Ap = 


sin y + 


A 


(8. 5.4) 


tan y + B 

where A and B are constants obtained from fitting to Eq. (8.5.2) for dry air. For 
A = 0.00095.B = 0, the error is reduced to 3 mm at 20 °, and 1.2 cm at 15°. 

The most accurate mapping for every path can be obtained by a ray trace 
program, but as long as low elevation angles ( Y < 15°) are not of interest, the 
above semi- empirical formula is adequate to the accuracy of about 0. 3 cm. 

It should be noted that the dry refractivity profile used in ray trace pro- 
grams is carefully determined from fitting to the radiosonde balloon data (see 
Glossary) [8. ?]. The fluctuations of the dry profile are small, and should not 
cause significant mapping errors for elevations higher than 15° [8.3], [8.4]. 

In conclusion, it should be mentioned that ionospheric calibration is not 
required when laser ranging hardware is utilized in the proposed multilatera- 
tion ranging schemes. 
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In conclusion, it should be mentioned that ionospheric calibration is not 
required when laser ranging hardware is utilized in the proposed multilateration 
ranging schemes. 

8.6 TRAJECTORY GENERATION USING GEOMETRIC PRINCIPLES 

During the operational phase of performing the type of multilateration 
experiments discussed in this report a rough knowledge of the orbital ephemeris 
(See Glossary) will be required. It is emphasized that the ephemeris of the 
orbit is not required for the determination of the station locations; it is only 
required for the purposes of predicting suitable viewing windows (Para- 
graph 8. 2) and for purposes of aiming the laser beams in the approximate direc- 
tion of the satellite. 

A rather unique side benefit of the multilateration scheme is that a very 
precise ephemeris of the satellite can be generated as long as the satellite is 
visible to at least three stations. Specifically, if the latitude, longitude and 
elevation of three stations are assumed to be known*, then the satellite ephem- 
eris can be determined fl. l] [4. 7] (See also Appendix D). The process of 
trilateration can be used to determine the position and velocity of the spacecraft 
to a high degree of accuracy over a short time span (5-10 minutes); conventional 
orbit determination methods would require a much longer duration to achieve the 
same accuracy (1-2 days). 

8. 6. 1 Difference Between Trilaterated and Dynamically Determined Orbit 

Usually, the objective of orbit determination is to obtain the fundamental 
constants of the orbit. Hence, a solution is sought for the six constants or 
elements, say, the position and velocity vectors at some specified time, which 
then can be used in conjunction with the second order differential equations of 

*This assumption obviously is made also when conventional orbit determination 
schemes are used. 
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motion and, upon integration, yield the future position and velocity of the 
satellite. To accomplish this, a model must be assumed; i.e., the potential 
field of the primary (See Glossary) must be specified, the perturbations due to 
Sun, Moon and planets must be included, and various other factors, e. g. , 
radiation pressure, drag perturbations, etc. need to be represented. 

On the other hand, when an orbit is tr ilaterated , the model is defined by 
only a set of equations which describe the intersection of spheres. (See Sec- 
tion 7). There are however a set of parameters common to both techniques 
which will be discussed presently. 

The difference between the orbit determination scheme and the trilatera- 
tion scheme can be discussed further by means of Figure 8-4. 



TRUE ORBITAL 
PATH (DOTS) 


Figure 8-4. Polar View of the Orbital Path of a Trilaterated 
Satellite. T, and a Dynamically Determined Satellite, D 
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This Figure illustrates the paths that would be predicted if the orbit determination 
method is trilateration, T or, dynamic determination, D. The dynamic method 
gives rise to a mean or smooth curve which is a best fit to the assumed model 
discussed above. The trilaterated path yields what is called an osculating 
curve, i. e. , a tangent curve as close as possible to the actual conditions exist- 
ing at that instant of time. It is obvious that since fewer error sources exist, 
the ephemeris generated by trilateration techniques will be closer representa- 
tion of the instantaneous positions and velocities at a given time. 

This statement is best explained with the aid of Figure 8-5. 


(a 3 , S 3 ) 



Figure 8-5. Schematic Orbit Determination Process 

Standard orbit determination procedures will use many observations 
a^, 6^ taken over a pass of the satellite, and by means of the assumed model 
reduce all of the observations to a set of initial elements at some adopted epoch. 
Schematically, if six observations a^, 6 , a^, 6 a y 6^ are obtained, then 
using standard techniques the observations which are separated in time by 
(t^ - t^), and (t^ - t^) can be mapped into the elements at t^, i. e. , the positions 
and velocities (r, i*) ^ can determined. But these elements are a function of 
the model for the potential field, perturbations, etc. , and are thus only defined 

For example, a , the right ascension, and, 6, the declination. 
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with respect to the assumed model. These elements are therefore mean 

elements in some sense which is not easily defined; however the error sources 

peculiar to the model are contained in the determined elements. When the 

orbit is trilaterated, many of these error sources and the process of 

obtaining a mean fit to the dynamic model are circumvented. Hence, when the 

T 

position and velocity (r, r) , where T denotes tr ilater ation, are obtained 

*2 

via tr ilateration the computed elements will certainly be closer to the true 
values of position and velocity. This is important because, as illustrated in 
Figure 8-6, when a prediction to a future state of the satellite is desired the 
initial conditions which will be assigned to the predictor equations are closer 
to the true values and therefore the predicted path should follow the true path 
more closely. 



Figure 8-6. Orbital paths as a Function of Initial Conditions 
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8. 6. 2 Listing of Error Sources 

Table 8-1 lists the error sources peculiar to each scheme and those 
sources common to both methods. 


Table 8-1 
Error Sources 


Dynamic Method 

Geometric 

Method 

Common 

• Gravitational Potential 

• Simultaneous 

• Measurement of Data 

Effects 

ranging 


• Sun, Moon, Planet 


• Station Coordinate 

Perturbations 


Locations 

• Drag/Lift 


• Constant of Precession/ 

Perturbations 


Nutation 

• Radiation Pressure 


• Sidereal Time 

Effects 


Determination 

• Magnetic Effects 


• Smoothing of Data 

• Other Small Effects 



(known and unknown) 




The error sources common to both methods are listed under the assump- 
tion that the trilaterated orbit will be transformed from the geometric coordi- 
nate system to the standard inertial coordinate system (See Section IV). For 
purposes of orbit evolution this is a necessity because an inertial coordinate 
system wherein the equation of motion are usually evaluated must be used for 
purposes of determining the future orbital position and velocity. If the equa- 
tions of motion are written in an Earth fixed, i. e, , rotating system, no advan- 
tage is gained because the Coriolis components would have to be added to the 
equations. 
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8, 6. 3 Outline of Tr ilater ation Technique 

From an operational point of view, after a relatively short duration, the 
following smoothed range polynomial (Paragraph 8. 3) will be available: 

p. = a. + b. (t - t R ) + c. (t - t R ) 2 , i= 1 , 2, 3, (8. 6. 1) 

1 1 1 X\ 1 JL\ 

where i=l, 2, 3 represents each station, and t^ is the maximum rise time for 

.K 

the three stations. A direct differentiation yields the range - rate as: 

p. = b. + 2 c i (t - t R ), i= 1 , 2 , 3 , 

so that six constants, i. e. , p^ and pu are available at time t (the time at which 
the satellite position and velocity are desired). 

From the known values of the station locations, namely, the latitude, 
longitude and elevation, and the other parameters displayed in the Common 
column of Table 8-1, the inertial position of the stations is determined by 
the methods of Appendix J. Then utilizing the standard trilateration technique 
[1.1] , the position and velocity components are determined in the inertial sys- 
tem defined in Section IV, i. e. , in symbolic terms; 

p i> P jy ^ jy ^ * x, y, z, x, y, z, t, 

where x, y, z, x, y, z are the satellite coordinates. The process can be 
repeated as often as desired over the time duration of common visibility to all 
three stations. This is accomplished by rejecting the data associated with the 
earliest strike of the satellite, accepting data from the latest strike and 
obtaining new values of a, b, c in Eq. (8. 6. 1) and (8. 6. 2) via the process of 
least squares. 
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Figure 8-7. Sequential Processing of Data 

Figure 8-7 is provided to clarify this technique. Suppose that data 
has been fit via least squares over an interval of duration A. As the illustration 
shows, nine data points would then be used to yield the values a, b, c of the 
curve fitting procedure. By means of tr ilateration it follows that the psuedo 
observations (p, p ) at strike 5 can be obtained from the fit and subsequently 
(r., r)j_ can be computed. Suppose that the satellite state is desired at strike 6. 
Since the curve fit interval or duration will be restricted, the best procedure 
to follow is to reject the data at strike one, pick up the data at strike 10, 
i. e. , use an interval of duration B = A, and refit (process sequentially). The 
computation with a, b, c now corresponding to the data in interval B will now 
yield (r, r)^ ; etc. By this method the curve fit polynomial will be of an 
osculating nature and thus provide the maximum accuracy. 

If a prediction of the satellite position and velocity at a time of non- 
visibility is desired then the position and velocity at some adopted epoch can be 
used along with the equations of motion and by numerical integration the desired 
state can be determined. 
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8. 6. 4 Relative Data Accuracies 


It is shown in Appendix F that the standard deviation of the data, i. e. , 
range, can be obtained as: 

3 

% = 2 /N V < 8 - 6 - 3 ) 

where N is the total number of data points and Up^ is the standard deviation of 
the actual data from the assumed model (curve fit). It is a simple extension 
to show that [l.l]: 


2 

cti = — — tr > 

p y/N T pf 


(8.6.3) 


where cr^ is the standard deviation of the range-rate as obtained via differen- 
tiation and T is the total duration over which measurements are taken. 

With these errors, a standard mapping can be performed and the errors 
in the satellite positions and velocities can be obtained. 

8. 6. 5 Conclusion 

The discussion of this Section demonstrates that a highly accurate method 
of predicting the future state of a satellite is a natural by-product of the multi- 
lateration technique. In brief, this makes the direct dependence on large soft- 
ware orbit determination programs a cir cumventable problem. 
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WEATHER CONSTRAINTS 


8. 7 

One of the important constraints limiting the reliability with which a 
station can successfully range to a satellite during a given pass is weather. 
Fortunately, an area of exceptional interest to geophysicists involved in 
earthquake hazard estimation is the American Southwest, including northern 
Mexico. In this region the weather is often clear, and obeys predictable 
patterns. There, the strategy of laser satellite geodesy suggests that the 
establishment of benchmarks be avoided in a given area during its season of 
unfavorable weather. Bad weather is discussed herein under the following 
headings: Nimbus (raincloud), and Stratus and Cumulus clouds, associated with 
the two types of rainy season in the Southwest. Then cirrus (icecloud); fog, 
both as mist and as icefog; and dust are discussed. 

8.7. 1 Nimbus, Stratus, and Cumulus 

The rainy seasons in the American Southwest are determined by the 
presence or absence of a tropical continental (cT) air mass which forms in 
summer over northern interior Mexico and adjacent parts of the U. S. , and 
which dissipates in winter. It appears on weather maps as an oval area of 
very low pressure with the axis running typically NW-SE through Yuma, 

Arizona. In summer and fall, when this air mass is present, Pacific storms 
are diverted to near the U. S. -Canadian border. In this season, tropical mari- 
time air (mT) from the Gulf of Mexico brings storms to New Mexico, to Arizona, 
and occasionally as far as Salt Lake City. In winter and spring, when this air 
mass disappears, some storms from the stable polar maritime air mass of 
the north Pacific pass over the American deserts, bringing rain across the 
California coast and to the Mojave, Chihuahuan, and Magdalena Deserts as 
far south as La Paz, Mexico. Thus, this stable summer low pressure area 
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divides the American Southwest broadly into three regions. In California west 
of the Sierra Front, and in Baja California north of 30° N, rain and clouds are 
frequent from December to April. Such cloud formations persist for several 
days and nights. In Arizona, New Mexico, and northern interior Mexico, brief 
thunderstorms are frequent from June through September. The storms themselves 
last several hours at most, but they are accompanied by cumulus clouds by day 
and various kinds of broken clouds by night which persist over long periods. 
Between these two regions, in a long corridor stretching from Death Valley 
across Baja California below 30° N to the Pacific Ocean, rainfall averages 
less than 12 cm. per year. Actually, the averages are misleading because 
some years have no rainfall, while in wet years the precipitation comes for 
the most part in sudden cloudbursts. In this third region, there is no rainy 
season. 

8. 7. 2 Cirrus 

High, feathery cirrus would impair laser ranging observations even 
in dry weather. They are characteristic of the passing of storm'fr ont s , and 
therefore are most frequent in the rainy seasons described above. 

8.7.3 Fog 

The west coastal deserts of Baja California are visited frequently by 
fog during spring and summer, as far south as 28°N. Fog is frequent on the 
California coast during fall and winter. A stable high pressure area often 
forms over eastern California and Nevada in winter, bringing fog to the San 
Joaquin Valley and icefog to central Nevada which persists for days, even in 
nominally clear weather. 
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8.7.4 Dust 


Storms of sand and dust can be expected in the American Southwest 
wherever there is loose, dry soil. April is the windiest month and therefore 
the dustiest. Winter brings comparative atmospheric calm and freedom from 
dust storms throughout the Southwest. 

8. 7. 5 A Model of System Reliability 

Atmospheric transparency, then, is not typically limited by extensive 
weather systems, such as are so common east of the Rocky Mountains, which 
cloud out the whole Southwest. Instead, it is probable that an event which 
closes down one station will leave the others unaffected. In this case, a fair 
model of the reliability of the laser ranging system is given by the statistics 
of the binomial distribution. Given n stations, each of which has probability 
p of acquiring useful ranges on a given pass of a satellite'', the probability that 
exactly k stations will acquire useful ranges is 

B (k, n;p) = k , ( n J_ k) , p k (l-p) n " k . (8.7.1) 

Then the probability that at least k stations will obtain useful ranges is 
obtained by summing B (k, n;p) from k to n: 


P (k, n;p) 



n! j , , , n-j 

P J (1-P) • 

j! (n-j)! 


(8. 7. 2) 


The function P (k, n;p) is tabulated in Tables 8-2 and 8-3 for the 6- and 
the 4-station cases, respectively. One obvious advantage of the 4-station 
method (using intercontinental stations; Section I) is that the probability of 

5,: Strictly speaking each station should be assigned a different probability 
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success for the whole system is higher, given stations each of a certain 
reliability. The most cost-effective strategy, however, may be to use 8 or 9 
stations, giving high probability that every satellite pass will give sufficient 
data for a 6-station solution. 

In conclusion, weather offers no serious problem to laser satellite 
geodesy, at least in the American Southwest. Further study would be needed 
to state system constraints for other areas. 

Table 8-2 

The probability, given n stations, each of which has probability p of 


acquiring useful ranges on a given pass of the satellite, that at least 6 stations 
will acquire useful ranges. 


\ 

P 

. 750 

. 800 

. 850 

. 900 

. 950 

n 






6 

. 178 

. 262 

. 377 

. 531 

. 735 

7 

. 445 

. 577 

. 717 

. 850 

. 956 

8 

. 679 

. 797 

. 895 

. 962 

• 994 

9 

. 834 

. 914 

. 966 

.992 

. 999 

10 

. 922 

. 967 

. 990 

. 998 

1. 000 
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Table 8-3 


The probability, given n stations each of which has probability p of 


acquiring useful ranges on a given pass of the satellite, that at least 4 stations 
will acquire useful ranges. 


\p 

. 750 

. 800 

. 850 

. 900 

. 950 

n 






4 

. 316 

. 410 

. 522 

. 656 

. 815 

5 

. 633 

. 737 

. 935 

.919 

. 977 

6 

. 831 

. 901 

. 953 

. 984 

. 998 

7 

. 929 

. 967 

. 988 

. 997 

o 

o 

o 

r— H 

•II 

8 

. 973 

- 990 

.997 

=1. 000 

=1. 000 
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8.8 SAFETY CONSIDERATIONS 

The use of powerful laser beams for the multilateration system 
contemplated herein presents some minor problems with respect to possible 
accidental injury. Although the energy levels preclude the possibility of damage 
to objects inadvertently crossing the laser beam, such as: birds, helicopters, 
airplanes, etc. , there is still the distinct possibility of injury to the human 
eye; specifically, to the retina. 

According to Laser Focus [8.5] a panel of the American National 
Standards Institute has come up with the following tentative figure for the maxi- 
mum exposure of laser light permissible that does no harm to the eye. Specifi- 

_ 7 

cally, for visible light the maximum permissible exposure (MPE) is 5x10 

, 2 -9 - 5 

Joule/cm for an exposure time of 10 to 2 x 10 seconds. This means, 

_ 3 

that for an anticipated 0. 5 Joule laser beam spread into a cone of 10 rad it is 
only safe to look into the beam at the minimum distance of about 6 km. The 
previous statement can be verified by the following simple relationships: 


MPE 


5x10' 


Laser pulse ener; 
beam area 


0. 5 Joule 
tt 0 r 


where r is the critical distance from the laser and 9=10 rad is the beam 


divergence. Obviously special precautions for low flying airplanes have to 
be taken. 


Accordingly the most important and rather obvious safety measures 
are listed below. The implementation of these safety measures is necessary 
but not difficult to achieve. They are as follows: 

• Operation of the laser devices should be in a controlled area. 

• Operations should be restricted to authorized personnel. 

• A ban on spectators in the controlled area without "appropriate 
supervisory approval" and protective measures, should be imposed. 
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Accurate control of the beam's path must be achieved. 

"Special care" must be taken in the use of such viewing aids as 
telescopes, which may increase the hazard. An interlock or filter 
must be provided. 

Eye -protection devices "should be used as supplements to 
engineering controls whenever the beam path is not completely 
enclosed." Devices should be "specifically designed for protection 
against radiation from the laser system being used. " 
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SECTION IX 


LASER RANGING HARDWARE STUDY 



The requirements on a satellite ranging system that would be suitable for 


use with the overall concepts presented in this report are stringent, i.e., it 
must determine the slant range to an Earth orbiting satellite to an accuracy of 
at least ±1 cm. Present laser ranging systems have uncertainties at least a 
factor of 10 larger than this. There is, however, a general concensus that 
the accuracy stated above is presently attainable [9. l]. 

This section presents a discussion of laser ranging systems, and the 
requirements necessary to upgrade their performance in order to attain ±3 cm 
or better accuracy. 


9. 1 PRESENT SYSTEMS 

9. 1. 1 Present Approach to Precision Laser Ranging 

Techniques and apparatus have been developed to determine the range to 
an object by measuring of the time- of- flight of a pulse of an optical wavelength 
[9.2]. The availability of high brightness pulsed lasers, capable of pro- 
ducing short duration optical pulses, has provided the basis for high preci- 
sion ranging systems. These systems closely resemble microwave radar in 
their mode of operation, i.e,, they operate on the principle of taking a simple 
measurement of the time required for a pulse to travel to and from the trans- 
mitter and a selected target. Laser transmitters have three distinct advantages 
over microwave ranging devices, namely, higher peak power, shorter wave- 
length, and shorter pulse duration. 

Many time-of-flight laser ranging systems have been constructed, mostly 
for field use by the military. The application being considered here, i.e., 
ranging to near Earth satellites, has not seen any commercial manufacturing 
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activity as of this date. Several government research centers such as 
Goddard Space Flight Center and Smithsonian Astrophy sical Observatory and 
the Air Force Cambridge Research Laboratory have constructed laser systems 
for ranging to near Earth satellites [9. 3], [9.4]. The University of Maryland, 
under contract to NASA, [9. 5] has assembled a system for ranging to the Moon. 
Research organizations in France and Russia also have satellite and lunar 
tracking stations using lasers. 

A general description of these present laser ranging systems will be 
given, along with a summary of their performance parameters. This will be 
followed by a discussion of the sources of errors and finally by some comments 
on the possibility of reducing these errors. 

9. 1. 2 Present Laser Satellite Ranging Systems 

The laser satellite ranging systems currently in existence all operate 
in basically the same manner. A general description will be presented first, 
followed by a summary of significant differences in the systems. Frequent 
reference will be made to Figure 9-1, which shows a block diagram of a laser 
satellite ranging system. Only the salient features are presented; the descrip^ 
tion of ancillary equipment used to provide automatic operation is omitted. 

The laser used for all existing systems discussed in Paragraph 9. 1. 1 is 
of the ruby type, and has been designed to make the device emit short intense 
pulses of light. These Q-spoiled, or giant pulse, ruby lasers typically produce 
pulses that are 10-30 ns in duration. Each pulse contains a few Joules of 
energy. After the laser fires, a small portion of the transmitted pulse is 
picked off by a detector which produces a pulse to start the timing electronics. 
The outgoing laser pulse goes through appropriate optics to give it the desired 
beam divergence angle, and to aim it in the proper direction. The object 
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Figure 9-1 Block Diagram of a Typical Laser Satellite Ranging System Current in U 






being tracked usually carries optical cube-corner retroreflectors in order to 
enhance the intensity of the return signal. The return pulse is collected by the 
receiving telescope and is focussed through a filter upon a photomultiple r tube 
(PMT) (see Glossary). The electrical pulse from the PMT furnishes the stop 
pulse to the timing electronics. A range gate is included in this signal path, 
in order to prevent unwanted signals from prematurely stopping the timer. 

The measured time delay is then read out, and the system is reset for the next 
event. 

The performance parameters of several existing laser tracking systems 
are summarized in Table 9-1. Two of the stations, GSFC and SAO, are used 
to track Earth satellites for geodetic purposes. The other two are used 
primarily for tracking the retroreflectors left on the Moon by the Apollo 
astronauts. 


9. 1.3 Sources of Errors in Present Systems 

Present satellite laser ranging systems are limited in the accuracy they 
can attain by several factors. The major cause of uncertainty is produced by 
the duration of the transmitted pulse and the resolution of the timing electronics. 
Minor uncertainties are the result of uncalibrated delays in the station elec- 
tronics and in the atmosphere. The laser ranging systems mentioned in 
Table 1 generally have overall accuracies that are in the 15-30 cm range. 

Clearly, the minimum resolution of the timing system will limit the 
overall accuracy. With a resolution of 1 ns, the range uncertainty would be 
±15 cm on a single measurement. This error can be reduced by performing 
repeated measurements and averaging the results. This technique is compli- 
cated by the rapid rate-of-change of the satellite's range, particularly for low 
Earth orbiting satellites. 
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Table 9-1. Performance Parameters of Several Laser Ranging Systems 
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The duration of the transmitted pulse can produce timing errors on the 
order of the pulse duration if the amplitude of the received pulse is small 
(a few photons). However, errors due to this effect will be relatively small 
(=1 ns) if the return pulse amplitude is large enough to permit "leading edge" 
detection in receiver electronics. 

Uncalibrated delays in the station electronics can produce errors on the 
order of one nanosecond. The major contribution to this error is the transit 
time delay variations in the receiver photomultiplier tube. As mentioned 
above, the start and stop pulses for the timing electronics are derived from 
two different photodetectors. The start pulse typically comes from an 
avalanche photodiode or phototransistor while the stop pulse comes from a 
PMT. The transit time of the PMT usually is 40 ns, and must remain constant 
during the course of a tracking pass. Small temperature changes or PMT 
power supply changes can produce changes in this delay, and hence cause a 
change in the internal station delay. 

The extant cube corner bearing satellites are all equipped with arrays 
of retro- reflector s . The physical extent of these arrays produces a smearing 
of the return signal, leading to errors of the order of the size of the array. 

In addition, the array of reflectors causes extensive amplitude fluctuations 
in the return signal due to interference effects (see Appendix I). This can lead 
to range errors even if the receiver electronics does compensate for the signal 
level fluctuations, since the return pulse shape may be severely distorted. 

The Earth's atmosphere produces an additional delay of its own, in 
addition to refracting the laser beam for non-zero zenith angles. The effects 
of this delay and refraction are for the most part well understood (see 
Paragraph 8. 5). Small variations exist due to the normal thermal turbulence 
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prevalent in the atmosphere. These variations will lead to errors generally 
significantly less than a nanosecond in magnitude. 

9. 2 IMPROVED APPROACH TO LASER RANGING 

The utility of a satellite laser ranging system which has significantly 
smaller overall error performance characteristics has been clearly estab- 
lished in the earlier Sections of this report. There is a real need for laser 

g 

ranging systems with accuracies approaching a few parts in 10 . This cor- 
responds to a timing accuracy of 0. 2 ns for typical Earth satellite ranges. 

The question can certainly be asked; Is it possible, using present technology, 
to produce a laser ranging system which is nearly a factor of 10 better than 
existing systems? The answer, as JPL envisions the present systems, is 
"yes". The components required to upgrade the system accuracy presently 
exist, and such a system can be constructed. 

The sources of error previously mentioned, i. e., pulse duration, timing 
resolution, and uncalibrated station delays can all be reduced significantly, 
thus leading to improvements in ranging accuracy. Two recent developments 
have led to the availability of equipment which will aid significantly in reducing 
the first two error sources; these are the mode locked laser and the high resolu- 
tion interpolating time interval meter. The third error source can be virtually 
eliminated by arranging the optics so that both start and stop pulses originate 
in the same PMT. These improvements are discussed further in the following 
Subsections. 

9.2.1 Laser 

The most significant advance toward eliminating a source of error is the 
development of the "mode locked", Q-spoiled laser. These devices are capable 
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of producing a single pulse of light of duration shorter than 0. 1 ns, which has 
a total energy in the pulse of more than 0. 5 Joule (5 gigawatt peak power). 

These lasers can be either of the ruby (0. 694 pm) or neodymium 
(1. 06 pm) type, with the appropriate Q- spoiling modifications to effect a giant 
pulse mode of operation. The addition of a mode locking device, usually a 
saturable absorber which modulates the laser cavity loss at its fundamental 
longitudinal mode frequency, produces a phase locking of all the oscillating 
longitudinal modes. This results in the production of several (usually 5-10) 
very narrow pulses during a period of 30 ns. It is necessary to include an 
electro-optic shutter outside the laser cavity in order to provide a means 
for releasing only one of the several mode locked pulses. 

Such a mode locked laser can produce pulses of 20 millijoules energy. 
The addition of two or three stages of amplification raises this energy to 
0. 5 joule or more. There are some differences between the ruby and 
neodymium lasers, the most notable being that the neodymium units are more 
powerful. A summary of performance parameters of one commercial laser is 
presented in Table 9-2. 

9.2.2 Timing Electronics 

Devices to measure the time interval between two pulses have existed for 
a long time. The typical unit counts the number of clock pulses which occur 
between the start/stop pulse pair in question. Of course, the clock pulse 
generator is well stabilized to minimize the effects of clock frequency errors. 
The timing error is at least ± one clock period. Such units typically have a 
10 MHz clock, which leads to a ±100 ns timing error corresponding to a range 
error of 15 m. 
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Table 9 _ 2. Mode Locked Ruby Laser Characteristics 


Manufacturer 

KORAD Department 
Union Carbide Corporation 
2520 Colorado Avenue 
Santa Monica, CA 90406 

Laser Type 

Q- switched, mode locked Ruby 
plus three-stage amplifier 

Pulse Energy 

0. 6 joules/pulse 

Pulse Width 

100 - 200 picoseconds 

Repetition Rate 

10 pulses per minute 

Beam Divergence 

2 milliradians full angle, 1/2 energy 

Cost 

$55, 000 

The single stage oscillator is Q- switched by a flowing dye cell, 
which also produces the mode locking action. A single pulse from 
the series of mode locked pulses is extracted using an electro-optic 
shutter, and amplified by a three-stage amplifier. The oscillator 
and amplifiers are operated from two separate power supplies, and 
cooled with a closed cycle heat exchanger. 


Recently, timing devices have been built which have much higher 
resolution than previous units. They still operate on the same basic idea, 
counting pulses with the aid of a relatively low frequency clock, but they are 
able to determine the time of occurrence of the start and stop pulses by inter- 
polation between clock pulses. 

Three different commercial manufacturers are currently producing time 
interval counters which operate on this principle. All of these have minimum 
resolutions of about 0. 1 ns. A more detailed investigation of these counters 
was undertaken; the results are reported in Paragraph 9. 3. 2. 
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9. 2. 3 Timing Pulse Pickoff 

A description of how timing errors produced by the manner in which 
the start and stop pulses for the timing electronics are derived has been 
presented previously. Since the goal of this section is to outline a system 
that is capable of timing accuracies of 0. 1 ns, it will be necessary to discuss 
a method for minimizing the errors caused by the drifting of station internal 
delays. 

By far the best approach toward solving this problem is to arrange for 
the same photomultiplier to serve as the detector for both the start and stop 
pulses [9. 6]. The procedure would be to arrange for a small amount of the 
transmitted pulse, only a few hundred photons, to reach the receiving PMT. 
This will produce a pulse at the output of the PMT which is used to initiate the 
timing sequence. The optical signal from the satellite received by the PMT 
will produce another pulse at the PMT output, which will serve to stop the 
timer. Since both the start and stop pulses have to traverse the same electron 
multiplier chain, cables, etc., only a few milliseconds apart, the effect of 
variations in the internal delays of all these devices will be minimized. 

The use of the receiver PMT as the source of the start pulse is beset 
with two complications: electromagnetic interference (EMI) from the laser, 
and backscattered signal from the transmitted pulse. The EMI produced by the 
firing of the laser, electro-optic shutter, and amplifiers may cause a false 
or premature signal at the PMT output. An additional signal may be produced 
by atmospheric backscatter and stray, off-axis radiation reaching the PMT 
during the first few nanoseconds after the transmitted pulse leaves the system. 
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Control of the EMI problem may be difficult. Suitable shielding, careful 
attention to eliminate ground loops, and application of microwave absorbing 
material such as Eccosorb will be required. The backscatter can be eliminated 
by the use of a mechanical shutter in front of the PMT. The shutter would open 
one millisecond after the transmitter fires, and remain open until after the 
return pulse is received. 

9.2.4 New Design for a Retroreflector 

The extant retroreflector bearing satellites, a total of six, all have 
large arrays of small diameter retroreflectors. A short duration laser pulse 
will be smeared out in time when reflected from such an array, unless the array 
happens to be oriented perpendicular to the line of sight. For precision 
ranging with short duration pulses, a different retroreflector arrangement would 
be necessary in order to reduce the smearing effects. 

A single aperture cube-corner type retroreflector would be the best. 
However, the relatively narrow effective field-of-view of a cube-corner would 
prevent ranging from several stations simultaneously [9. 11]. A. special array 
of 4 or 8 cube-corners could be constructed, with the apexes of all cubes 
positioned so that they are nearly coincident with each other. This device 
would produce reflections from almost any angle of incidence. The individual 
cube-corners would be required to have a relatively large aperture, 5 cm 
diameter or more, and would have to be accurately constructed with the three 
faces mutually orthogonal to high accuracy. The return beam diameter at a 
receiving station site would be small enough so that the velocity aberration 
caused by the satellite's motion will shift the return beam completely off the 
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receiving telescope aperture. It will thus be necessary to somehow spread the 
return beam slightly. Further research would be required to develop a suitable 
cube-corner retroreflector array with the properties suitable for simultaneous 
use with several stations. 

9. 2. 5 Performance Potential 

The discussion presented in the previous sections has indicated that 
several of the major sources of errors in laser ranging systems may be 
eliminated or materially reduced. If a system were built utilizing these 
improvements, what would be the remaining sources of error, and what would 
be its expected absolute accuracy? With the major error sources eliminated, 
some of the lesser error sources can be discussed. 

There are several clearly defined sources of small errors, and 
probably many that we have not yet identified. Four of the more predominant 
effects are: 

• photomultiplier tube response time, 

• pulse height variation from shot-to-shot (see Appendix I) 

• calibration of atmospheric delay and refraction (see Section 8. 5) 

• atmospheric turbulence effects. 

The first two items listed are likely to be the more important of the 
lesser errors. They have been the subject of an extensive experimental 
investigation: the results are reported in Paragraph 9.3. 1. The general conclu- 
sion is that: commercially available photomultiplier tubes and pulse height 
compensation apparatus exist, but need to be improved. 

The problem of calibration of atmospheric delay has been extensively 
treated by several groups. The most detailed of these was performed by 
Hopfield [8. 1]. The conclusion reached was that; with good knowledge of the 
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atmospheric temperature and pressure at the laser station site, it is possible 
to estimate the atmospheric delay to within a few millimeters at the local 
zenith, and to within 2 centimeters at an elevation angle of 20° . Clearly, 
experimental verification of these conclusions is necessary, but the results 
are encouraging since the resultant errors are small (See Paragraph 8. 5 for a 
more detailed discussion of the problem of atmospheric calibration). 

Atmospheric turbulence effects on the overall system accuracy should be 
small, since their time average will be zero. However, scintillation and 
multipath effects may deteriorate the performance by producing confusing 
signals in the PMT, and by causing the PMT output pulses to vary in amplitude. 

The deleterious effects of all these lesser sources of errors will be 
small. However, in order to provide a simple system with the highest possible 
accuracy potential, it is necessary that these effects be understood. No attempt 
has been made to experimentally investigate the atmospheric effects during the 
course of this one year study. Any follow-on program should include a 
detailed study of atmospheric effects, in the hope of further improving system 
accuracy. 

9.3 EXPERIMENTAL PROGRAM 

Two major areas were investigated experimentally in this study, namely: 
the effect of photomultiplier response time on the ranging accuracy was studied, 
and the accuracy and resolution of the commercially available high resolution 
time interval meters was accessed. As part of this study a small demonstra- 
tion laser ranging system was constructed in order to verify all the features 
necessary to achieve high accuracy range determination. 
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9.3. 1 Photomultiplier Tube Response 

If an optical ranging system based upon time -of -flight techniques is to 
have fractional nanosecond accuracy, it is clear that the majority of the com- 
ponents in the timing portion of the system must have response characteristics 
that are consistent with such a time requirement. In the discussions of equip- 
ment it has been shown that two of the critical components, the laser and the 
time interval electronics, have potentially adequate response times. * The 
same cannot be said of the photomultiplier tube. 

A typical PMT has a response time on the order of a few nanoseconds, 
with the duration of the entire transient lasting tens of nanoseconds. The 
transit time required from the reception of an optical pulse at the PMT 
photocathode to the detection of an electrical signal out of the anode terminal 
(at the base of the PMT) is usually 20 - 40 ns. The nature of the trajectories 
of the electrons as they travel from the photocathode through the electron 
optics (if any) to the first dynode, the subsequent trajectories through the 
electron multiplier section, the collection efficiency of the anode, and the 
way that the anode is connected to the world outside the PMT envelope, all 
contribute to deterioration of the response time of the PMT. Of all these 
effects, the first has the most predominant effect, accounting for about 
80 percent of the time smearing [9. 6]. 

Photomultiplier tube manufacturers are presently becoming aware of the 
need for improved transient response of their products. Their efforts have 


*The phrase response time in this report is defined to be the time required to 
reach a preset trigger level. The response time is related to, but not the 
same as, the rise time. After an event has passed a trigger level, the shape 
itself, or how faithfully it reproduces the initiating pulse shape is of no 
basic importance. 
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extended in many different directions, from development of radically new 
types of PMT's*, to the addition of coaxial feedthroughs for the anode signal 
leads. ** 

An extensive amount of research into the transient response properties 
of PMT's has been done by various organizations in the nuclear instrumenta- 
tion field [9. 7]. Their experience has led to a choice of a photomultiplier tube 
which is for the most part conventional in its construction, except that the 
large photocathode-first dynode drift space has been greatly reduced. 
Furthermore, the photocathode material used in this tube: gallium arsenide, 
is by far the best choice for use with a ruby wavelength laser. 

The photomultiplier tube chosen is a "Quanticon" type C31034, manu- 
factured by RCA. Table 9-3 presents a summary of its salient characteristics. 
The tube has a quantum efficiency that is 12 percent at the ruby wavelength; 

this is the highest efficiency of any current tube at that wavelength. It has 

7 

an internal multiplication gain of 10 at 2000 V, and an exceptionally low 
dark current. Figure 9-2 is a photograph of a sampling scope display of the 
single photo output pulse shape. The rise time is approximately 3 ns, with 
a width at half maximum of 5 ns. This response is exceptionally fast for a 
PMT of this type, that is, of more or less conventional construction. 

The more exotic types of PMT's mentioned previously, the c rossed-field 
type and the channeltron type, both have some undesirable features in their 
performance characteristics. The crossed-field PMT has an excellent transient 
response (<0.5 ns) but it has a low multiplication gain. The channeltron PMT's 
transient response is no better than a conventional PMT. They both have 

*The Sylvania crossed field PMT and the Bendix channeltron are examples of 
this type. 

**Both Amperex and RCA tubes use this feature. 
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Table 9~3. Characteristics of RCA C31034 Photomultiplier Tube 


Photocathode 

Material 

— Gallium Arsenide 

Spectral response curve no. 

- 128 

Wavelength of maximum response 

— 830 nm 

Typical quantum efficiency 

550 nm 

- 14% 

694. 3 nm 

- 12% 

860 nm 

- 9% 

Electron Multiplier 

Dynode material 

— Beryllium-Copper 

Number of stages 

- 11 

Structure 

— In line, electrostatic focus 

Dark Current Pulse Summation 

1/8 to 16 photoelectrons at 20° C 

- 10, 000 


the more conventional types of photocathodes, with the correspondingly low 
efficiency at the 0.694 pm ruby laser wavelength. These factors make the 
utility of these tubes questionable. In addition, the experimental set up used 
to test the PMT response utilized an infrared diode laser whose output wave- 
length is beyond the cutoff of these tubes. 

The characteristics of the gallium arsenide photocathode are so superior 
that a program should be considered to incorporate it into the more exotic 
photomultipliers mentioned. Another more conventional PMT, the C31024 
developed by RCA, has potential as a fast response time tube. It has a 
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Figure 9-2. Sampling oscilloscope display of the average shape of the 
dark current pulse from an RCA C31034 photomultiplier tube. The 
abscissa is 2 ns/major division, and the ordinate is 10 mv/major 

division. 


conventional photocathode, but only five dynodes instead of the eleven in the 
C31034. Each of these dynodes is fabricated from gallium phosphide, and has 
a single stage gain of more than 10. The five stage GaP multiplier has as much 
gain as the eleven stage conventional multiplier. The characteristics of a 
photomultiplier tube with a gallium arsenide photocathode and gallium 
phosphide dynodes would probably make it extremely useful as a sensitive, 
fast response detector for laser ranging applications. 


JPL Technical Memorandum 33-605 


9-17 




An improved version of the RCA C31034 already exists, and is 
designated C31034A. This PMT has a slightly different photocathode, with 
enhanced quantum efficiency in the red part of the spectrum. 

9.3.2 High Resolution Time Interval Meters 

One major requirement for a centimeter accuracy ranging system is a 
device that will determine the time-of-flight of the laser pulse to an accuracy 
of 0. 1 nanoseconds. A conventional time interval meter, which measures time 
by counting clock pulses between the start and stop signals, would require a 
10 GHz clock frequency in order to have a resolution of 0. 1 ns. Such devices 
presently do not exist, although they may in the near future due to micro- 
integrated circuit technology. 

Electronic timers with resolutions of 0. 1 ns have been developed 
recently, and are available from three different commercial manufacturers. 
These devices all operate on the same principle: utilizing a basic clock fre- 
quency of only a few megahertz; the time between start and stop signals is 
determined by counting the clock pulses, with the added ability to determine the 
time of occurrence of the start and stop signals between the clock pulses. This 
is accomplished by the use of electronic interpolators that are able to effec- 
tively resolve the basic clock period into 0. 1 ns intervals. 

Figure 9-3 shows a block diagram of an interpolating time interval meter. 
The pulse counter commences counting the first clock pulse after reception of 
the start trigger signal, and ceases counting with the next clock pulse after 
reception of the stop trigger signal. The start (stop) interpolators measure 
the time between the reception of the start (stop) signal and the next clock 
pulse.. The special purpose processor performs the mathematical addition and 
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subtraction necessary to produce the time reading for the display and/or 
output. A typical time interval meter utilizes a clock pulse generator with a 
100 ns period (10 MHz), and interpolators capable of resolving the clock period 
into 1000 parts, each of 0. 1 ns duration. 

The accuracy of these devices is limited by the usual problem of clock 
generator frequency stability, interpolator linearity, and input pulse shapes. 
For purposes of the present program, the frequency stability does not 
represent a problem since the total time duration being measured will only be 
a few milliseconds. The interpolator nonlinearity may be a source of error, 
but as reported below the effects of linearity can be measured. Errors due to 
input pulse shapes can be minimized by providing fast risetime, constant 
shaped pulses to the input of the timing circuitry. 

Three commercial units are available that have the 0.1ns resolution 
required for this program. The salient features of these high resolution 
counters are summarized in Table 9-4. 

The procedure for testing the accuracy of these counters is based upon 
the following assumption: the major source of error will be in the start and 
stop interpolators, with systematic errors in measurement of long time 
intervals being attributed to a clock frequency error. A test procedure was 
developed to determine the variation in measured delay caused by interpolator 
errors. No attempt was made to measure clock frequency errors. * 

Figure 9-4 shows a diagram of the experimental apparatus used to 
generate the start and stop pulses that are used as inputs to the time -interval - 
meters under test. The pulse generator is set to deliver a so-called NIM 
Standard fast pulse. The pulse has a duration of 10 ns, an amplitude of 

*Clock frequency stability will not be a problem, since it is always possible to 
provide an external clock source with the required stability. 
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Table 9-4. High Resolution Time Interval Meters 



JPL, Technical Memorandum 33-605 
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-1 volt, and a risetime of 1 ns. It is split by a power dividing tee. One of the 
split pair is terminated in 50f2 and serves as the start pulse; the other is 
introduced into a delay line which is terminated in 50 ft and serves as the stop 
pulse. The pulse generator is set to run at a frequency which is hot harmoni- 
cally related to the internal clock of the time-interval-meter under test. This 
provides essentially a random time of occurrence of the start pulse with 
respect to the internal clock pulses. 

Repeated single measurements of the time delay between the pulses was 
made, and a record of each occurrence kept. The average value of the 
measurements will be a reasonable indication of the delay between the pulses, 
and the standard deviation will be an indication of the errors introduced by 
the interpolators. * 

The results of the measurements are summarized in Table 9-5. The 
standard deviation of the measurement of each of the three timers indicates 
that each one was performing within manufacturer specification. That is, all 
three units had a resolution of about 0. 1 ns, and the data are in agreement. 

If an interpolator was not performing in a linear fashion, the random time of 
occurrence of the timing pulses would produce a larger spread of the measured 
events, and hence a larger standard deviation. 

The conclusion of this test of commercial time -interval- meters is that 
they have adequate accuracy performance to satisfy our timing requirements. 
The present implementation calls for a unit that has a common input, due to 
the fact that the start and stop pulses originate within the same PMT, and will 
travel along the same transmission line. The Hewlett-Packard unit is the only 

*No corrections were made for the effects of transmission line loss on the 
shape of the stop pulse. This effect is a constant, and will not alter the con- 
clusions about interpolator errors. 
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Table 9~5. Measurement of Time-Interval Meter Interpolator Error 


Unit Tested 

HP 5360A 

EG & G 


Nanofast 536B 


+ 5379A 

TDC 100 


+ A- 3 plugin 

Measured 

28.6 - 

2 

17.25 - 

8 

5. 9 - 1 

Times and 
Number of 

28. 8 - 

58 

17. 375 - 

72 

6. 0 - 10 

Occurrences 

28. 9 - 

223 

17.50 - 

15 

6.1-69 


29.0 - 

196 

17.625 - 

5 

6. 2 - 20 


29. 1 - 

18 





29. 2 - 

2 





29.7 - 

1 




Number of 
Events 

500 


100 


100 

Mean Value, 

28. 936 


17.397 


6. 108 

ns 






Standard 

Deviation, 

0. 085 


0.079 


O'. 058 

ns 







one presently that has a common input. The other manufacturers indicate that 
they can provide such an option on their units. 

At the present time, our recommendation would be to use the Hewlett- 
Packard H12-5360A plus the HOI-5379 plugin. This recommendation is based 
solely on cost. The H12 option eliminates unneeded features, and is priced 
at $5445. 00. 

9. 4 DEMONSTRATION LASER RANGING SYSTEM 

An experimental apparatus was devised to demonstrate the response 
of a photomultiplier tube to short duration light pulses, and to demonstrate the 
concept of ranging with only a single photodetector. The apparatus utilizes a 
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pulsed injection diode laser as a source of light, a sampling oscilloscope to 
display the signals, and a high resolution time-interval-meter to measure the 
transit time of the light pulses. Figure 9 “5 is a photograph of the apparatus. 

9.4. 1 General Description 

Figure 9-6 presents a diagram of the apparatus, showing the general 
layout of the various optical components. The light from the diode laser is 
collimated into a 20 mm diameter beam by a lens. * This beam is then split 
into two beams of approximately equal intensity. One beam goes to a nearby 
mirror and the other travels to a cube-corner retro- reflector (CCRR) located 
some distance away. The light from each reflector returns to the beam splitter, 
half of each beam being directed toward the photomultiplier tube. A narrow 
band interference filter is placed in front of the PMT to cut down on the stray 
light that reaches the photocathode. An aperture of 8 mm diameter limits the 
PMT field of view. Photographs of the optical component arrangements are 
presented in Figure 9-7. 

The duration of the light pulse that the laser emits is very short, only 
3 ns. Thus, two distinct light pulses reach the PMT, with the time difference 
between them being related directly to the difference between the optical path 
lengths of the paths the two beams traverse. The nearby mirror produces the 
start pulse for time interval meter, and the distant CCRR produces the stop 


*The lens, a 30 mm focal length Hastings triplet, had no astigmatic correction 
to accommodate the effective shape of the diode laser's output aperture. 
Consequently, the beam had an approximately elliptical shape in the far 
field. 
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Photograph of the Demonstration Laser Ranging System, Showing the 
al Assembly, Laser Pulse Generator, Power Supplies, Constant 
Fraction Timing Discriminator, and Time Interval Meter. 






9-27 


JPL Technical Memorandum 33-605 






9-28 


JPL Technical Memorandum 33-605 


Figure 9-7 (a). Optical Assembly with Covers in Place 
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Figure 9-7(b). Optical Assembly with Covers Removed. The laser and collimating lens are in the 
box at the bottom. The large box at left-center contains the beam splitter, mirror, and start 
pulse attenuator. The box at right-center contains another attenuator and a narrow band 
interference filter. The photomultiplier tube is in the closed box on the extreme right. 



pulse. An optical attenuator* placed in front of the nearby mirror allows 
control over the intensity of the start pulse, and another optical attenuator 
placed in front of the PMT allows the intensity of both pulses to be adjusted 
simultaneously. The electrical signals from the PMT are fed to the constant 
fraction timing discriminator, a device that produces standard shaped output 
pulses with a constant internal delay that is independent of input pulse 
amplitude. The output of the constant fraction discriminator triggers the 
time interval meter. 

9.4.2 Laser Pulse Generator 

The arrangement for providing a short duration electrical pulse to the 
diode laser is shown in Figure 9-8. The Tektronix Model 109 pulse generator 
is a charge line type, with the duration of the pulse determined by the length 
of the 50 £7 charge line (pulse duration = 2x electrical line length), and the 
amplitude of the pulse controlled by the amplitude of the voltage stored on 
the charge line (pulse amplitude = 1/2 stored voltage). The repetition rate is 
approximately 250 pps. A signal pickoff is placed in the output line of the 
device, to permit monitoring of the laser driving pulse, and to provide a 
synchronizing signal for a sampling oscilloscope. (Fischer [9.8] describes 
a similar method for driving diode lasers.) 

The laser diode, an RCA TA7867 gallium aluminum arsenide type, emits 
light at a wavelength of 850 nm. This type of diode used has a threshold current 
of 10 amp, a maximum current of 25 amp, and a peak power output of 5 watts. 
The laser is mounted on a copper heat sink, with a 49. 8G resistor connected 
in series with the center conductor to provide proper termination for the 


*A11 attenuators used in this experiment are of the Schott glass neutral density 
filter type. 
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Figure 9-8. Pulse Generator Circuit for Providing Short Duration, Fast Rise 
Time Pulses to a Gallium- Aluminum- Arsenide Diode Laser. 




electrical drive pulse. A photograph of the laser in its mount is shown in 
Figure 9-9. 

This method for driving the laser has the distinct advantage of providing 
a fast risetime, short duration pulse to the laser. Figure 9-10 is a photograph 
of the sampling oscilloscope trace of the laser pulse. The pulse duration was 
3 ns, and its amplitude was set to just over threshold, usually 12 amps. The 
major disadvantage of this method is that the charge line voltage must be 
excessive. Most of the energy in the electrical pulse is dissipated in the 
termination resistor at the diode mount. The termination is necessary to pre- 
vent the signal at the diode terminals from reflecting back into the pulse gen- 
erator and to prevent the voltage across the diode terminals from reversing. 
Note that reverse polarity pulses will permanently damage the laser. 

To provide a 12 amp pulse to the laser, the signal generator must deliver 
a 600 volt pulse, which requires 1200 volts on the charge line. The Tektronix 
Model 109 pulse generator is not designed to run normally in this manner. 
However, its operation has been made quite reliable by doing two things; 
namely, providing a large resistor (2Mf2) in the charging path from the external 
power supply, and using a charge line on only one side of the generator's SPDT 
reed switch. 

9.4.3 Photomultiplier Tube Circuit 

The photomultiplier tube, an RCA C31034, is mounted and wired in a 
nearly conventional manner. Figure 9-11 presents the wiring diagram. 

Figure 9-12 is a close-up photograph of the socket wiring. The dynode voltages 
were provided in the usual manner, from a resistive voltage divider chain, as 
recommended for this tube by RCA. The last four dynodes were bypassed with 
appropriate capacitors, to prevent saturation from large pulses. 
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Figure 9-9. Gallium- Aluminum- Arsenide Diode Laser Mount. The Lase 
RCA TA786 7, is Mounted in a Cylindrical Copper Block. A 49.8 £2 
Terminating Resistor is Installed Inside the Center Hole of the 
Block, in Series With the Center Conductor of the Diode. 
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Figure 9-10. Sampling Oscilloscope Trace of the Electrical Pulse 
Driving the Diode Laser. The abscissa is 1 ns/major division, 
and the ordinate is 100 volts /major division. 


Extreme care was taken in the connection of the output signal cable to 
the PMT socket. Hard copper shielded 50£2 coax was used, with the shield 
ending 3 mm from the pins on the tube socket. The last two dynodes were con- 
nected to the coax shield with 100 pf high frequency capacitors*, in order to 
maintain some semblance of a transmission line configuration through the tube 
socket. (The C31034 is supplied with internal strip leads between the anode/ 
last dynodes assembly and the tube socket pins. This arrangement improves 


^Microwave ceramic capacitors with metal strip leads supplied by the American 
Technical Ceramics Co., Huntington Station, New York. 
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HV INPUT 


ELECTROSTATIC SHIELD 
AROUND TUBE ENVELOPE 



Figure 9-11. Photomultiplier Tube Socket Wiring Schematic. 
The tube is an RCA C31034. 
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Figure 9-12. Close-Up Photograph of the Photomultiplier Tube Socket Wiring, Showing 
Details of the Output Cable Connection to the Tube Socket Pins. 


the high frequency characteristics of the tube.) Electromagnetic interference 
from the pulse generator was minimized by placing a block of Eccosorb in the 
PMT housing. 

A time domain reflectometer was used to check the tube/socket assembly 
for mismatches on the signal output cable. The result is displayed in 
Figure 9-13. The strip line inside the PMT appears to have a 75 £2 character- 
istic impedance. 

Improvejnents could be made in the PMT wiring. RCA suggests that the 
dynode voltages can be modified to provide improvement in the saturation 
characteristics at large pulse amplitudes. A matching network could be 
placed between the tube socket and the 50£2 output cable, to provide a better 
match to the internal strip line, or the external cabling and circuitry could be 
changed to a 75Q characteristic impedance. No serious problems in the 
operation of this demonstration laser ranging system were encountered, there- 
fore no attempt was made to include these improvements into the system. 

9.4.4 Constant Fraction Timing Discriminator 

The constant fraction timing discriminator is basically an active pulse 
shaping network. It produces a constant shaped output pulse, independent of 
the shape of the input. The discriminator is intended for use with scintillation 
and other particle detectors in timing systems for nuclear instrumentation, 
but is ideally suited for this present application. Its main function is to prevent 
the stray, low amplitude PMT output pulses caused by background and dark 
current from triggering the timing electronics. The discriminator must do so 
in such a manner that it does not introduce any undesirable delay variations 
of its own. The device has the built in ability to correct for input pulse height 
variations, greatly reducing the "walk" error. 
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The constant fraction timing discriminator used for this series of 
experiments was the ORTEC model 463. Its salient characteristics are 
summarized in Table 9~6. 

Table 9-6. Constant Fraction Discriminator Characteristics 


Manufacturer: 

ORTEC, Inc. 

100 Midland Road 

Oak Ridge, Tenn. 37830 

Model No. : 

463 — Constant Fraction Discriminator 

Input : 

Accepts negative pulses to 10 V without 
saturation 

Dynamic Range: 

Walk error less than ±150 psec for a 
range of 100 mv to 10 v 

Output; 

NIM Standard fast negative output 
signal 

Cost: 

$550. 00 
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9.4.5 Experimental Results 


The demonstration laser ranging system was used to study the performance 
of the photomultiplier tube and timing system with signal levels similar to those 
expected in an operational laser ranging system. The response of the photo- 
multiplier tube to pulses of different amplitude was studied, and the absolute 
accuracy and random errors of a complete ranging system were investigated. 

9. 4. 5.1 Photomultiplier Tube Response 

The shape and time of arrival fluctuations of the electrical pulses out of 
the RCA C31034 PMT were studied by recording the pulse shapes on a sampling 
oscilloscope as a function of the intensity of the optical pulse reaching the 
photocathode. Figure 9-14 is a photograph of a typical sampling oscilloscope 
trace of signal from an RCA C31034 PMT. The particular PMT used, serial 
number Y07716, produces a 10 millivolt output pulse from one photoelectron 
when the operating potential is 2000 volts. The trace in Figure 9-14 is 1.4 volts 
at the peak and has a duration that is similar to the single photoelectron 
response shown in Figure 9-2. The assumption is made that the peak of the 

ajc 

transient scales linearily as the number of photoelectrons in the signal. The 
signal was produced by 140 photoelectrons, or assuming 10% quantum efficiency, 
1400 photons. This number checks with the calibration data of the laser and the 
filters used to attenuate the signal. 

Figure 9-14 can also be used to estimate the transit time dispersion. 

Since a sampling oscilloscope was used to obtain the photograph, the resultant 


*The justification for this assumption is based upon the fact that the duration of 
the laser pulse, 3 ns, is less than the impulse response of the PMT, including 
transit time dispersion. Thus almost all the photoelectrons and their 
secondaries reach the PMT anode in a time that is less than the single 
photoelectron impulse response time. 
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Figure 9-14 Sampling oscilloscope trace of the electrical signal from an 
RCA C31034 PMT. The optical pulse from the laser was 
3 ns long, and approximately 1400 photons reached the 
photocathode. The abscissa is 2 ns /major division, and 
the ordinate is 0. 5 v/major division. 
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trace is effectively an average of many pulses, 140 per major division. The 
time spread in the leading edge of the pulse at the half peak value is approxi- 
mately 0.2 ns. The transit time dispersion has a Poisson distribution, and if 
the 0.2 ns is assumed to be a 3cr limit, then the standard deviation of the transit 
time is approximately 0. 068 ns. 

Similar data was taken for several different signal levels, with the result 
plotted in Figure 9-15. The standard deviation of the response time dispersion 
is 0. 1 ns at a signal level of 100 photoelectrons. This is an acceptably small 
response time variation, and will not contribute significant errors to a range 
measurement. The minimum acceptable signal should produce at least 100 photo 
electrons at the PMT photocathode. 

9. 4. 5. 2 Experimental Ranging System Accuracy 

The entire experimental apparatus, as shown in Figure 9-6, was used in 
conjunction with a Hewlett-Packard (HP) 5360A/5379A high resolution time 
interval meter to determine the overall system accuracy. The optical signal 
levels were adjusted to correspond to those expected in an operational system, 
with 100 photoelectrons /pulse being the standard acceptable signal level. 

The HP 5360A/5379A time interval meter was used in conjunction with 
another HP accessory, the HP 5375A keyboard, which gives the time interval 
meter the capability to perform some computations from an internally stored 
program. The unit was programmed to calculate the mean value and standard 
deviation of a successive series of time interval measurements from the rang- 
ing experiment. The number of samples in the ensemble could be set to any of 
the following numbers: 1, 10, 100, 1000, or 10,000. 

The program presents the result in distance, rather than time. The value 

Q 

for the speed of light that is used, 2. 997141 x 10 m/sec, is for dry air at 
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NUMBER OF PHOTOELECTRONS/PULSE 


Figure 


9-15. Variation of the Response Time of an RCA C31034 
PMT with Signal Level. The PMT was Operated 
with a Divider Network Shown in Figure 9. 11, and 
a Supply of 2000 Volts. 
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mean temperature and pressure that exists at the test range used for these 
measurements. The correction for changes in atmospheric conditions is 
negligible because of the short ranges used, typically less than 50 meters. 

The measured values of mean range and standard deviation are presented 
in Table 9-7, for various values of range and for different signal levels. The 
range was adjusted by placing the cube-corner- retroreflector at three different 
distances from the apparatus. Each distance was carefully measured with a tape 
measure. The signal levels were adjusted by inserting neutral density filters 
in front of the PMT and/or in front of the nearby mirror. The amplitude of the 
start and stop pulses could thereby be adjusted independently. 

The results indicate that the system was capable of measuring the range 
to within an accuracy of 20 mm, provided the start and stop signals were the 
same size and contain at least 100 photoelectrons. The standard deviation of 
the measurements indicates that random errors in the system are no larger than 
expected. The data for the longer ranges shows an increase in the standard 
deviation above the 20 mm level. This is due to severe atmospheric turbulence 
effects, caused by the optical path being too close to the ground surface, and was 
further complicated by the small diameter receiving aperture. A larger diam- 
eter receiving aperture would materially reduce this effect. It is pointed out 
in Paragraph 8. 5 that atmospheric turbulence effects for reasonable elevation 
angles are of the order of 10 mm or less. 

There is an error which occurs when the start and stop pulses are not the 
same size. This error has been traced to the ORTEC Constant Fraction Dis- 
criminator in our experimental setup. Subsequent tests of it indicate that it 
does in fact make such errors, which were reduced further by performing the 
so-called "walk" adjustment. Further research will be required with constant 
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Table 9-7 Experimental Ranging System Accuracy 
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Return signal levels variable due to atmospheric turbulence effects. 










fraction discriminators in order to verify that such devices can function 
properly without introducing appreciable errors. 

9.5 Transportable Laser Tracking Station 

The basis for the entire concept of precision determination of interstation 
distances and coordinates is a transportable laser station capable of accurately 
measuring the range to a cube corner bearing Earth satellite. This section 
discusses the requirements that a suitable station must meet, and describes 
a station in moderate detail. A cost estimate for development and procurement 
of a network of these stations is presented. 

9.5.1 System Range Equation 

Utilizing the standard communication range equation [9.4], with appro- 
priate corrections for the laser beam angle [9. 1 0] , the useful energy at the 
photomultiplier photocathode is defined by Eq. (9.5. 1), as 


N 
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2 2 
\ Dj 
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2 2 4 

hce t £ e t R 


^a^r 1 ^ 


q P t 


(9.5. 1) 


where the various terms and their units are presented in Table 9-8, and suit- 
able values for these various terms are also listed. Substitution of the typical 
values into Eq. (9.5. 1) yields 

N = 814 photoelectrons. 

This is a sizable signal. Note that one of the design criteria is that the return 
pulse should contain at least 100 photoelectrons. The system could operate at 
a range of 1.7 Mm and still have a suitable return. 
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Table 9-8 Description of Terms and Typical Values for 

Equation (9.5.1) 




Description and Units 

Typical Value 

Number of photoelectrons liberated from the 


PMT photocathods 


Receiving aperture diameter (meters) 

0. 5 

Retroreflector diameter (meters) 

0. 05 

Slant range (meters) 

10 6 

Laser operating wavelength (meters) 

0. 694 x 10" 6 

Product of Plank's constant and the speed of 

1. 986 x 10" 25 

light (meter-joules) 


Transmitted beam divergence angle at the 

10‘ 3 

half power points (radians) 

\Tz 

Retroreflector return beam divergence 

io -4 

angle at the half power points (radians) 

n/2 

Transmitter optics efficiency 

0. 8 

Atmospheric transmission efficiency 

0. 7 

Retroreflector efficiency 

0. 5 

Receiving optics efficiency (including 

0. 25 

filter) 


Detector quantum efficiency 

0. 10 

Energy of transmitted pulse (Joules) 

0. 5 
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It is thus clearly feasible to obtain reasonable signal returns with a 
modest system (1/2 joule pulses, 50 cm collecting aperture). 

9. 5. 2 Tracking Mount Description 

The basic specifications for this assembly are that it be capable of 
tracking an Earth orbiting satellite in a 500 km or higher orbit, with full sky 
coverage to 20° elevation angle and an accuracy of ±0. 5 milliradian. The 
optical package includes 50 cm diameter receiving optics and -a 4 cm diameter 
transmitting optics. Special requirements include the need to be fully auto- 
mated for unmanned operations and the need to be transportable to remote 
stations with a minimum setup time and installation crew. 

The ground station that best meets the above requirements as well as 
the criteria of minimum cost is described as follows. An X-Y mount will be 
placed on air ride shock absorbers in a specially built trailer. The trailer 
will also contain all of the laser equipment and control electronics required 
to operate in remote locations. The basic concept is that special pads would 
be prelocated and aligned at remote sites; the trailer would then be brought to 
the location and lowered onto the special pads, thus requiring only minimum 
alignment prior to operation. The basic power source would be commercial 
power or independent generators. 

The tracking mount will be an X-Y configuration with the X axis aligned 
with the trailer. Several types of mount configurations were studied, AZ-EL 
and polar mounts have gimbal-lock problems to get the required sky coverage 
and fork type X-Y mount has no advantages since bearing errors are not a 
severe problem. 
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The basic structural element for this mount will be standard 20 inch thin 


wall, seamless steel pipe. Major assemblies will be welded and stress 
relieved. Axis bearings will be ABEC class 5 which are available at a reason- 
able cost. The major assemblies will be insulated to minimize the thermal 
deflection errors. 

The axis of the mount will be driven directly with DC torque motors of 

g 

50 lb-ft capacity (6. 8 x 10 dyne-cm). Axis rates are low because of the X-Y 
configuration, only l°/sec velocity and 0. 02 °/sec acceleration. The servo 
electronics will consist of a control chassis which will be integrated with the 
station computer and a drive power amplifier. Current, rate, and position 
feedback loops will be closed to insure the stability of the mount and to minimize 
the tracking error. 

To provide digital readout of axis position an encoding system will be 
geared to each axis. Although some angular errors are introduced by this 
technique, there is considerable cost savings over a direct shaft encoding 
system. The encoder type selected will be geared 36 to 1 with the axis, the 
resolution will be 0.001° and the accuracy equal to 0.005°. Gearing errors 
will be less than 0.003° /axis. 

Receiving optics will utilize a 50 cm diameter f-2 reflector with a detec- 
tor package mounted at the prime focus. A flip mirror will be utilized to 
provide a Newtonian visual system for alignment purposes. The mirror 
materials will be standard, there is no apparent advantage to use replica or 
aluminum mirrors. The laser assembly will be stationary off the X-axis and 
a series of four mirrors will be provided to transmit the output beam up to the 
Y-axis. Field alignment equipment will be provided to verify the parallelism 
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of the output beam and the receiver optics. Figure 9-16 shows a diagram of the 
optical paths in the mount. The start pulse path would be provided by a fibre- 
optic bundle, with appropriate attenuation to control the amplitude of the start 
pulse. Figure 9-17 is an artist's sketch of the optical mount assembly. 

The pointing accuracy of the tracking mount is calculated by taking the 
root of the sum of squared peak values of the independent error sources. 

Table 9-9 presents the significant error sources. The 3<r value of the peak 
errors listed is 42 arc sec. On the basis of engineering analysis of the proposed 
mount design the principal error sources will be: axis orthogonality, mount 

alignment and the encoders. These errors could be reduced by increasing the 
cost of the mount, but it is felt that the present design is a satisfactory balance 
between cost and performance. Some reduction can also be achieved by calibra- 
tion of systematic errors and programming corrective offsets into the pointing 
computer. This programming complexity is not costed and the error calculation 
in Table 9-9 assumes no compensation. 

9. 5. 3 Return Pulse Detection System 

The detector assembly mounted at the focus of the 50 cm collecting optics 
must incorporate a shutter, narrowband interference filter, photomultiplier 
tube, and pulse discriminating electronics. 

The shutter is provided to prevent scatter light from the transmitted 
pulse from reaching the PMT. It must be a fast acting, electromechanical 
unit which opens a millisecond or so after the transmitter has fired. 

The narrowband interference filter is provided to permit day time opera- 
tion of the system. The bandwidth should be less than 1 nm, and it should have 
a peak transmission of at least 50% (such filters are readily available). 
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Table 9-9. Tracking Mount Errors 


Error Source 

Type 

Peak Value 
(Arc Sec) 

1. Telescope Tube def. 

Systematic 

0. 9 

2. Y-Axis bending 

Systematic 

3. 0 

3. Y-Axis bearing def. 

Systematic 

1 . 0 

4. X-Axis Torsion 

Systematic 

0. 9 

5. Y-Axis Bearing Runout 

Random 

7. 7 

6. X-Axis Bearing Runout 

Random 

3. 3 

7. Orthogonality 

Systematic 

20. 0 

8. Temp. Gradient in Telescope Tube 

Random 

6. 0 

9. Temp. Gradient in Base Legs 

Random 

2. 0 

10. Y-Axis Encoders 

Systematic 

18. 0 

11. Y-Axis Encoder Gearing 

Random 

10. 0 

12. X-Axis Encoders 

Systematic 

18. 0 

13. X-Axis Encoder Gearing 

Random 

10. 0 

14. Mount Alignment 

Random 

20. 0 


The photomultiplier tube should be one of the new Quanticon type (see 
Paragraph 9. 3. 1), in order to utilize their enhanced sensitivity at the ruby laser 
wavelength. Special care must be taken with the design of the PMT circuit, 
in order to provide an optimum match between the output of the tube and the 
external electronics. The circuit design also must prevent gain saturation in 
the PMT for large amplitude pulses. 

The last component of the detection system is the pulse height discrim- 
inator. This device is used primarily to discriminate between PMT output 
pulses that are caused by the return signal, and pulses caused by dark current 
and background radiation. The unit's output is a standard NIM fast negative 
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pulse, and occurs only when the input pulse has exceeded a preset height, 
presumably determined by the background and dark current levels. 

As discussed in Paragraph 9.4.4, this unit must introduce no variable 
delay of its own. The so-called constant fraction pulse height discriminator is 
designed to prevent the most common such variable delay, that caused by 
different amplitude signals. Although the constant fraction discriminator used 
in the experimental part of this program did not perform entirely satisfactorily, 
it is felt that such units can be constructed. 

9. 5. 4 Laser 

The type of laser to be used in this ranging system is clearly very 
necessary to the operation of the system. The laser must emit optical pulses 
that are of very short duration: 100-200 picoseconds. The mode-locked, 

Q-spoiled ruby laser discussed in Paragraph 9. 2. 1 fulfills this requirement, and 
has enough energy in each pulse to permit the use of modest size cube-corner 
retroreflectors and receiving apertures. 

At the present time, the KORAD Department of Union Carbide Corporation 
is the only laser manufacturer that claims to have such lasers "off the shelf". 
There are many other manufacturers of pulsed ruby lasers, and these should 
be contacted if there is an interest in pursuing this multilate ration concept. 

9. 5.5 Laser Ranging System Transporter 

A special purpose trailer is included as part of the ground station. 

Although transportability was not a prime requirement of this preliminary 
design, it was felt that the most cost effective solution to the various problems 
of transport, on-site set-up, electronic equipment housing, and protective 
covers could be solved with this special trailer. An artist's sketch of the 
assembly is shown in Figure 9-18. 
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Figure 9-18 Artist's sketch of a complete transportable laser ranging system, showing 
the optical mount in its transporter. The transporter is positioned over a concrete 
pad and the mount is lowered onto alignment pins. The transporter contains all 
the necessary station electronics, but does require an external power source. 


The mount and laser assembly will be carried in the trailer on air ride 
shock absorbers to minimize the shock of highway transport. On arrival at a 
ground station the trailer will be positioned over a specially designed concrete 
foundation. The foundation will have 3 pre- aligned register pins which mate 
with the base of the X-Y mount. The mount is lowered onto the register pin by 
deflating the air ride shocks. After the mount is set on the pins only a small 
amount of alignment verification tests will be required. 

The trailer will be equipped with a roll-back type of cover which is 
removed when the mount is in use. This type of cover is not as effective a 
thermal shield as a telescope type astrodome, but is a much more practical 
solution to this problem. 

The trailer also contains the space for the electronics and the laser power 
supplies. Also provided is the equipment air conditioning for the four racks of 
equipment. 

9. 5.6 Laser Ranging System Block Diagram 

A block diagram of the ground station is included as Figure 9-19. A small 
general purpose computer is used to run the station; including angle pointing, 
laser control, station monitor, data storage, range gating and many other 
r equi r ement s . 

It is planned that with an absolute minimum of information, such as the 
orbit parameters and laser firing times, the station could complete a mission 
automatically. The minimal input/output data is necessary because of the 
assumed poor quality of land communications to remote sites. 

The final sizing and configuration of the computer and its peripheral 
equipment will depend on the monitor and control requirements, which in turn 
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are dependent on the requirements and desirability of unmanned operations. 

For the purpose of this estimate a small computer with 1 OK of memory was 
used as a baseline. 

The frequency and timing equipment is dependent on where the stations 
will be used. The network requirements are that the four station clocks of a 
network be synchronized within 10 microseconds and the network be synchro- 
nized to Universal Time within a few milliseconds. 

Two clock synchronization systems are available. The first is the best 
and the least expensive, it depends on being within range of a commercial 
television station. Currently the NBS is broadcasting on the west cost for 
4 hours /day; a 1 megahertz signal and a time code as part of an experimental 
program called the TV line- 10 system. This system will provide a clock 
synchronization of 1 to 2 microseconds as well as a digital time code. By 
using this system crystal oscillators can be used at the sites. 

If the sites are not within range of a TV station then a rubidium frequency 
standard should be provided at each site. Also a master clock will have to be 
circulated among the sites on a weekly basis to maintain the 10 microseconds 
synchronization; the drift rate of a rubidium being about 1 microsecond/day. 

As part of the station equipment a remote weather station is included 
to monitor wind velocity and direction, ambient temperature, atmospheric 
pressure, and humidity. These parameters will be digitized and stored in the 
computer for remote readout. 

9.5.7 Cost Summary 

A cost summary for construction of a six station network of transportable 
laser ranging stations was prepared. The costs are broken down into three 
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areas: design, station equipment, and network support. The design cost 

presumes that commercial drawing and documentation practices are adequate, 
rather than the more elaborate JPL/NASA standards. The mount, laser, and 
equipment costs make no assumptions about quantity discounts. It is felt that 
the cost estimates are reasonable and possibly conservative. 

A network of six stations could be designed and constructed for a total 
cost of $1.85 million. This figure does not include site preparation costs. No 
attempt has been made to estimate the yearly cost associated with the operation 
of this network. This figure is very subjective at this time, depending upon 
the total number of sites necessary for geophysical significance, the ease of 
remote automatic operation, etc. 

Table 9 -10 presents the design cost estimate, broken down into the man- 
months required. Table 9-11 summarizes the cost of equipment and assembly 
for one complete station. Table 9-12 is an estimate of the equipment costs 
associated with support of the six station network. 


Table 9-10 Design Cost 


Item 

Engineering 

Drafting 

Mount 

3 mm 

6 mm 

Optic s 

1 mm 


Servo 

4 mm 

4 mm 

Encoders 

1 mm 

1 mm 

Trailer 

1 mm 

1 mm 

F oundation 

1 mm 

2 mm 

Alignment 

2 mm 


Programming 

3 mm 


T otal Man months 

16 

14 

T otal Cost 

$ 1 00K 
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Table 9-11 Equipment Cost 


Item 

Cost 

Mount 

$5 OK 

Optics 

$15K 

Servo 

$13K 

Encoders 

$ 7K 

T railer 

$25K 

F oundations 

$ 1 OK 

Assembly 

$25K 

Freq and Time 

$20K 

Computer 

$20K 

RCVR etc. 

$ 1 OK 

Laser 

$55K 


$25 OK 

Contingency 10% 

$25K 


$275K /Station 

T otal for 6 stations 

$1. 65M 


Table 9-12 Network Support Equipment 


Item 

Cost 

Transportable Primary 

25K 

Time Standard 


Alignment Equip. 

10K 

Test Equip, and 

15K 

Special Tools 


Spare Parts 

5 OK 

T otal 

100K 
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SECTION X 


RECOMMENDATIONS AND FUTURE ACTIVITIES 
| In the final Section of this report, four areas are discussed; i. e. , the 
recommended multilate ration system is outlined, a tentative future plan of opera 
tions is presented, a detailed but brief listing of the future areas of investigation 
is listed, and some pertinent recommendations are listed. No detailed attempt 
to plan the future of this project is presented herein; only sufficient information 
to focus on future activity is detailed. It is felt that the final program plan devel 
opment is of sufficient magnitude to warrant analysis and is, therefore, more 
properly included in the advanced systems planning phase (See Paragraph 10. 2). 
10. 1 RECOMMENDED SYSTEM CONFIGURATION 

Analysis performed within this report has treated the multilateration 
problem in a general context. The techniques developed herein can be applied 
to an arbitrary number of stations (24) and trajectory points. In this section, 
the practical "minimum" system configurations are examined. 

As has been demonstrated in Section VII, in a range-only system, the 
4- station/6-trajectory point configuration is a theoretical minimum with respect 
to the number of stations. On the other hand, no matter how many stations are 
used, at least 4 trajectory points are still required, with the 6-station/ 
4-trajectory point configuration being the minimum system in this class. 

Evidently, the cost of stations is an important consideration in designing 
a prototype system, whereas, the cost of additional measurements is not. 
Therefore, the 4-station system would have been the preferred configuration 
except for the very important fact that the 4-station problem is degenerate when 
the stations all lie on a plane (the coplanar degeneracy). As a result of this 
degeneracy, the 4-station system yields very poor sensitivities when stations 
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are nearly coplanar; this is inevitable if the stations are to be placed on the 
surface of the Earth over relatively short distances (1000-km or less). 


As discussed in Paragraph 7.4 and Appendix B the coplanar degeneracy 
persists even with 5 stations. Only when 6 stations are used can this degeneracy 
be removed. Hence, the analysis has dictated the following recommendations. 

• Primary Recommendation: ^-6-Station System 

This is a general-purpose system which can be made to work with 
short as well as long baselines. No satellite ephemeris informa- 
tion, beyond that required for aiming the laser beam, need be 
employed. Excellent sensitivities can be obtained for favorably 
located stations even if the stations are nearly coplanar. Hence, a 
6- station configuration is recommended as the foundation of a 
workable system . 

• Secondary Recommendation:' ■ — 4-Station System 

This system works only if stations are separated by large inter- 
continental distances, with one station located out of the plane 
defined by the other three by a large distance. Satellite ephemeris 
information is not required, but the altitude of the satellite must be 
high enough; i. e. , the altitude should be equal (approximately) to 
the Earth's radius. This will ensure that a sizable portion of the 
trajectory will be simultaneously visible from all four interconti- 
nentally separated stations. In turn this requires a high-power 
laser, and ranging at large zenith angles. Good sensitivities are 
possible if the above requirements are satisfied. 

Hence, if the cost of two additional stations is prohibitive, the four- 
station system is recommended. 
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Third Recommendation: 


3-Station, Collinear 


System 

This is a baseline -determination system utilizing three laser ranging 
stations approximately aligned in a straightline configuration. Due 
to the alignment requirement, only short baselines, less than 
200 km, can be determined. An approximate knowledge of the sat- 
ellite ephemeris (with an accuracy of about a meter) is needed to 
calibrate the misalignment error (see Appendix C). Accurate solu- 
tions are obtainable for at least one of the baselines. 

This is a minimum-cost system and is recommended (See Para- 
graph 10. 2) as an interim step in the overall project plan. Certainly, 
if this system is implemented, the construction of an additional 
station would be a trivial expense. The fourth station would then 
permit the 4-Station System to be exercised, etc. 

10.2 PROGRAM PLAN 

A schedule for the development of the proposed laser ranging system is 
shown on the next page. The schedule indicates the estimated time period (in 
quarters of a year, starting from the initiation of funding) during which various 
tasks will be performed. 

The first phase of the program will be devoted to advanced investigations. 
During the first six months, an intensive system analysis will be performed. 

This analysis will make considerable use of computer simulation, and will resolve 
all the questions raised in Paragraph 10. 3 except those specifically related 
to hardware. Laboratory tests to evaluate the critical hardware components, 
e. g. , lasers, PM tubes and time-interval meters, will be performed concurrently 
with the systems analysis. Certain geophysical questions (such as the required 
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Project Plan 



PRELIMINARY INVESTIGATIONS 

o) Advanced system analysis 

b) Final experimental evaluation of hardware 

c) Geophysical analysis and preliminary 
determination of suitable station locations 


SYSTEM DESIGN 

a) Hardware design 

b) Software development 


FABRICATION AND INITIAL TESTING OF 
PROTOTYPE STATION 

a) Fabrication of prototype station 

b) Testing of prototype station via horizontal 
ranging 


FINAL TESTING OF PROTOTYPE STATION 
AND POSSIBLE MODIFICATION 

a) Testing of prototype station via ranging 
to extant satellites 

b) Data reduction and evaluation of prototype 
station 

c) Possible modification of prototype station 
and possible deployment of a new, orbiting 
ref le c t i on -me chan ism 


PREPARATION FOR SYSTEM OPERATION 

a) Detailed geophysical analysis and final 
determination of suitable station locations 

b) Fabrication of 2 additional stations 

c) Testing 

d) Fabrication of 3 additional stations 

e) Testing of all 6 stations via horizontal 
ranging (with possible supporting tests 
via ranging to satellites) 


SYSTEM OPERATION 

o) Deployment of 6 stations in geophysically 
significant locations 

b) Extraction of raw data 

c) Data processing, system evaluation and 
geophysical interpretation 


10-4 


JPL Technical Memorandum 33-605 





frequency and accuracy of data measurements, and the selection of suitable 
station locations ' ) will be examined in detail in a sustained, moderate-level 
effo rt. 

System design will commence in the third quarter and extend through the 
fifth quarter. Detailed specifications for the hardware complement of a ground 
station will be developed. Considerable effort will be expended in order to 
achieve a hardware design which utilizes commercially available components 
and is of maximum cost-efficiency. A complete software package for the sys- 
tem will be evolved, debugged and finalized. Program efficiency will be 
emphasized due to the fact that an extremely large amount of data will have to 
be processed during the operational life of the system. 

Fabrication of a prototype station will begin during the latter part of the 
system-design phase. Any overlap between the fabrication and system design 
phases will be primarily devoted to procuring those long-lead items already 
predicated by the system-design study. Fabrication will be completed twelve 
quarters after the consummation of the system design. The prototype station 
will then be intensively tested by horizontally ranging to a ground-based 
retroreflector. Although these tests will suffer from the deleterious atmospheric 
effects common to all horizontal- ranging systems, they will provide a con- 
venient and inexpensive means for identifying hardware imperfections and esti- 
mating overall performance. 


* Initial investigations of suitable station locations have already been made, and 
are described in Section VI. 
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Final testing of the prototype station will commence in the seventh 
quarter. The prototype station will be used to track the extant, retroflecting 
satellites; these tests will not only provide a considerable amount of extremely- 
significant performance data, but will also discover possible operational prob- 
lems (such as difficulties in achieving proper satellite tracking). The data 
obtained from these tests will then be analyzed in detail so as to determine the 
statistics of the received- signal power, the severity of atmospheric effects, 
and -- most importantly -- the fluctuations in measured range due to different 
pulses being reflected from different corner reflectors on the satellite. An 
evaluation of system performance will then be made, and, if necessary, 
appropriate modifications to the system will be specified and implemented. If 
acceptable system performance is precluded because of the multiplicity of 
corner reflectors on the satellites, a new orbiting reflection-mechanism will 
be designed and launched. 

Once satisfactory results have been achieved from the prototype station, 
project effort will be focused on developing a complete operating system. An 

interesting system test can be performed as soon as three of the stations have 
been built, i. e. , in the twelfth quarter a test of the three station system 
will be undertaken (possibly with an airplane). 

A detailed geophysical analysis to select optimal station locations and the 
construction of three additional ground stations (identical to the prototype station) 
will dominate the final quarters of this project. All six ground stations will be 
tested via horizontal ranging to Earth-based retroreflector s ; also, if further 
evaluation of station performance appears necessary, each station will be tested 
via ranging to satellites. 
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The six stations will be deployed in their specified locations during the 
14th quarter, and the complete system will be operated for a period of a year. 
The considerable amount of data resulting from this operation will be processed 
and examined for geological significance, and a new program plan will be 
developed to schedule possible modifications and future deployments of the 
total system. 

10.3 FUTURE STUDY AREAS 

Paragraph 10. 2 indicates that the future program plan will commence 
with an advanced systems study. This study will be necessary to answer some 
of the details not considered under the present study phase and also to provide 
a second indenture analysis of some of the topics which were considered only 
superficially herein. Certainly, the main objective of this study, i.e., to 
prove the feasibility of the geometric technique, has been accomplished. How- 
ever prior to implementing the multilate ration scheme some further analysis 
will be required. These areas of investigation are briefly outlined in this 
Section as study tasks. From the length of the list it is obvious that a sizeable 
effort will be required to codify the remaining details. 

Task A-Systems Planning 

This will be a planning phase to develop a detailed milestone schedule 
and to set priorities on the remaining tasks which are defined herein. 

Task B-Solution Degeneracies 

A detailed study of degeneracies which occur when certain station and 
orbital configurations are assumed will be undertaken. This task which is a 
logical continuation of Appendix B of the present report, will identify poten- 
tially unstable configurations and explicitly define operational rules to circum- 
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vent these problems. It should be noted that these degeneracies are not fun- 
damental in the prevention of system operation; however, the attainment of 
maximum system accuracy demands a more detailed understanding of these 
effects. 

Task C-Number of Stations and Optimal Deployment 

This report has answered the question relative to the number of 
stations which are theoretically required in order to ensure nominal system 
operation. The advantages of using more than the minimum number of stations 
are not however understood. There is a possibility that use of more than the 
minimum number of stations will result in even more dramatic accuracy 
improvements than are displayed in Section V. This task will perform the 
analysis required to answer the question of what additional coordinate accura- 
cies are achievable via the use of more stations. The selection of stations for 
performing multilateration experiments will also be analyzed on a firm geo- 
physical foundation with much more detail than in Section VI. Once the theoret- 
ical configurations yielding optimum accuracy station locations are understood, 
a proper compromise with the geodetic constraints will be undertaken and an 
optimum realistic system configuration will be defined. 

Task D-Error Analysis 

Due to pressing project schedules the detailed error mapping of 
sources affecting station coordinate accuracies was undertaken using Monte 
Carlo techniques. In this task the explicit covariance mapping of the random 
error sources, and the effect of system bias errors will be performed. 

Task E-Use of Range-Rate Data 

The use of range -rate data was not considered in this study. A 
possibility exists that use of this type of data will result in further accuracy 
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gains. Task 4 would perform this type of analysis and incorporate any result- 
ing advantages into the overall system. 

Task F-Closed Form Solution 

There is theoretical justification to believe that the multilate ration 
equations can be solved in closed form. A solution to these equations in 
closed form would result in further understanding of the geometric process 
and in possible system simplification during real time operations. Analysis 
will be performed in this Task to verify the possibility of obtaining a closed 
solution to the fundamental equations of the process. 

Task G-Curve Fit Accuracies 

The detailed analysis to determine the maximum duration over which 
data can be fit to the adopted models will be undertaken in this task. This is 
a prerequisite to any operational system. 

Task H-Data Displays 

Useful displays of the parameters obtained from the multilate ration 
process will be defined. These displays will aid in the efficient transfer of 
information and in simplifying system operation. 

Task I- Pulse Shape Degradation 

A detailed study and verification of pulse shape degradation due to 
multiple reflections will be performed in this Task. 

Task J- Atmosphe ric Calibration 

An in-depth analysis of the best manner in which atmospheric calibra- 
tion can be performed in order to minimize measurement errors will be per- 
formed (see Paragraph 8. 5). This task is central to lowering the total system 
error budget. 

Task K-Tracking Analysis 
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Task K will deal primarily with defining the accuracy which will be 
required in satellite state estimation for purposes of pointing the laser assem- 
blies. As a function of these accuracies the suitable azimuth /elevation slav- 
ing equations will be derived and the respective hardware accuracies will be 
specified. 

Task L- Launch Vehicle Analysis 

The total retr oreflector system weight budget will be specified in 
this task and a systematic search of possible satellite /launch vehicles which 
can accommodate the retr oreflector package (as an additional mission objective) 
will be specified. 

Task M- Airplane Vehicle Analysis 

For purposes of interim experiments, the requirements imposed for 
implementation of Principle 3 (3 station system. Section I) upon suitable air- 
craft will be outlined. This will include hardware constraints and operational 
flight path maneuvers to ensure proper system operation. 

Task N-Three Station Calibration 

A detailed analysis of the effect of misalignment in the 3 station 
system will be performed. This error source was briefly analyzed in the 
present study. A better understanding of these effects is required. 

Task O-Inverse Trilateration 

The use of the method of inverse trilateration (Appendix D) for 
determination of new station locations as a function of previously defined 
station locations will be investigated for possible singularities and the best 
accuracies attainable in the solutions. The proper recommendations for 
operational implementation of this technique will be specified. 

Task P-Generalized Coordinate Systems 
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An investigation of the use of more generalized coordinate systems 
(See Section IV) will be performed. This task will be of fundamental importance 
if comparisons of parameters estimated via geometric and dynamic principles 
are to be undertaken. 

Task Q-Comparison of Geometric and Dynamic Techniques 

This task will concern itself with the information that can be 
extracted when residuals in the parameters obtained by means of geometric 
principles and dynamic principles are computed. There is theoretical justifi- 
cation which shows that certain astrodynamical constants can be calibrated 
using such methods. These constants would be identified and appropriate 
conclusions as to their accuracy would be drawn. 

Task R- Comparison With Other Systems 

The final task of the advanced study phase will concern itself with a 
detailed comparison of the final multilate ration system configuration with 
other competitive systems. Items such as accuracy, cost, ease of implementa- 
tion, etc. will be tabulated (See Section III). 

Task S- Software Development 

Throughout the duration of the advanced study phase, the development 
of the software tools required to analyze the measured data will be undertaken. 

Task T- Documentation 

The results of the advanced study phase will be documented in a con- 
tinuing effort such that at completion of this phase, detailed and final docu- 
mentation of the multilate ration technique will be available. This documentation 
will include the operational description of the non-real time software developed 
during this period. 
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10.4 PROJECT RECOMMENDATIONS 


In light of the previous discussion, it is recommended herein that the 
study phase of this project be continued on a moderate level; i. e. , a funding 
level of approximately $80, 000 to $100, 000 should be expended during the next 
study period for continued systems analysis. An equal amount should be utilized 
to answer hardware-oriented problems. About $50, 000 should be utilized to 
integrate the hardware and software subsystems, to provide future planning, 
and to develop appropriate system test plans. 

In summary, an effort with funding equal to approximately $250, 000 
would ensure continued effort on a system whose ultimate accuracy in the 
location of station coordinates can be measured in subcentimeter units. 
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SECTION XII 


GLOSSARY 


Bouguer Correction : The gravity values at an elevation h above sea level must 

be reduced to values at sea level. A number of corrections are therefore 
necessary. The correction for the effect of the mass in the layer between h 
and sea level suggested by Bouguer equals -2rrGph and is called the Bouguer 
Correction. Here p is the mass density and G the gravitational constant. 
Constant Fraction Timing Discriminator : An electronic device to measure the 
arrival time of a pulse accurately by integrating over its length, and taking a 
preset fraction of its height, so that the leading edge of the pulse can be 
determined. 

Corner Reflector: An arrangement of mirrors in such a fashion that any light 
beam will be reflected back independent of the attitude of the mirror system. 
Degeneracies : The process wherein the multilateration equations cannot be 

solved for station coordinates due to geometric or numerical instability in 
the equations. 

Drag : A force opposite to the velocity vector induced on a satellite or airplane 
due to air viscosity and density. 

Dynode : An electron current amplifier in a photomultiplier tube. Electrons 

are accelerated in an electric field, impinge on a plate (the dynode), and 

liberate secondary electrons. The gain of a dynode is between 2-3 and there 

are up to 14 dynodes in a photomultiplier tube. This makes the overall gain 
14 14 

equal to 2 -3 ; a rather respectable amplification. 
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Eccentric Anomaly: An auxiliary angle defined in terms of the true anomaly 


via the relationships: 


cos E 


cos v + e 
1 + e cos v 


sin E 


JT- 
1 + 


2 

e sin v 


e cos v 


where E is the eccentric anomaly, v is the true anomaly and e is the orbital 
eccentricity. 

Ecliptic : The plane of the apparent motion of the Earth about the Sun. 

Ephemeris: The position and velocity of a celestial body tabulated as a function 

of some independent variable, say, time. 

Epoch: A fixed time origin, i. e. , a starting time. 

Equinox : Usually called the equator and equinox, i. e. , the position of the 

Earth's equatorial plane and the position of some reference axis (x axis) at 
a given time. 

F ault: The surface of fracture between two blocks of rock which have moved 

relative to one another. 

Geodetic Datum : A surface, line, or point, and/or the measurements per- 
taining thereto, from which geodetic information may be calculated. 

Geodolite : The trade name of a commercial laser ranging device used for 

surveying manufactured by Spectra-Physics. 

Geometric Flattening of Earth : A measure of the difference between the polar 
and equatorial radius of the Earth (f - 1/298. 3). 

Gravimeter : An instrument (usually portable) for measuring the acceleration 
due to gravity. 
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Laser: Light Amplification by Stimulated Emission of Radiation. A devise to 

generate extremely monochromatic radiation making use of the coherence 
properties of stimulated emission. 

Loran System : An electronic navigational system by which hyperbolic lines of 
position are determined by measuring the difference in the time of reception 
of synchronized pulse signals. 

Magma Vent: The route by which magma, i. e. , molten rock under pressure 

is extruded through the Earth's crust, e. g. , into the neck of a volcano. 

Magnetic Lineations : A zebra-stripe pattern of the directions of magnetization 
in the rock of an ocean floor, generally found to be parallel to a mid-ocean 
ridge. The existence of this pattern is evidence both of changes in the direction 
of the Earth's magnetic field and of sea-floor spreading. 

Mode Locked : The internal modulation of a laser cavity leads to a radiation 

field which consists of many modes which are "locked" into oscillations with 
the same phase. 

Nutation: The perturbative effect upon the pole of the Earth caused by the 

Moon. 

Obliquity : Usually called the obliquity of the ecliptic; an angle between the 
planes of the Equator and Ecliptic (t = 23.° 5). 

Perigee : Point of closest approach of a satellite to the Earth; when referred 

to any other planet the corresponding term would be apoapsis. 

Plate Tectonics : Theory of: The hypothesis that the uppermost layer of the 

Earth's crust consists of a mosaic of rigid plates, and that oceanic trenches, 
continental drift, and at least some types of orogeny (mountain building) can 
be explained in terms of the relative motion of these plates. 
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Photodiode: A semiconductor devise, i. e. , a p-n junction. Minority carriers 
created by absorption of light within a diffusion length of the junction may dif- 
fuse to the junction and give a photocurrent. 

Photomultiplier : An extremely sensitive photon counting device which operates 
on the principle of the photo electric effect (release of an electron by a photon 
impinging on a solid state surface). Subsequent current amplification is 
accomplished by means of the release of secondary electrons. 

Primary : The mass about which local motion is occurring, e. g. , the Earth 

is the Secondary of the Sun (which is Primary) or the Moon is the Secondary 
of the Earth (which is the Primary); etc. 

Pseudo Observation : The observation obtained from raw data after performing 

mathematical /statistical operations on the raw data. 

Q -Spoiled: The Fabry-Pero cavity of a laser ideally has a very high Q. With 

the help of an electro-optic shutter or passive absorbing shutter, the Q will be 
considerably lowered to allow a light pulse to escape. The magnitude of Q is 
a measure of the sharpness of a resonance. 

Quartz Oscillator: A very stable electronic oscillator using a quartz crystal 
to define and stabilize its frequency via piezoelectricity. 

Radiation Pressure : Force per unit area exerted by electromagnetic radiation 

on a surface of a body. Its magnitude is (4it) ^ |E x H| which amounts to 

-5 -2 

4. 46 x 10 dyne cm for the Sun’s radiation at the distance of the Earth. 
Radiosonde : A balloon with attached measurement equipment to measure the 
temperature and pressure of the troposphere as a function of altitude. The 
information is telemetered back to ground by radio. 
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Range Gate: When a number of pulses are transmitted and subsequently 


received the possibility of confusion arises, i. e. , a pulse which is transmitted 
may not correspond to a pulse which is received. The range gate is elec- 
tronically present to the approximate arrival time and therefore rejects all 
other pulses and allows for an accurate measurement of the time of transit. 
Range Residual : Usually the difference between a pseudo observation and the 

raw data, i. e. , the difference between a nominal range and the raw range. 

2 2 

Range Square Differences : p^ - p^ where p^ and p^ are slant ranges evaluated 
at different points (times). 

Retroreflector : See corner reflector. 

Strain : The process whereby a given material undergoes a change in deformation. 
Stres s : The forces on a given area which give rise to internal tension or 

compres sion. 

True Anomaly : The angle, measured in the orbital plane of motion, between 
a line from the center of mass to perigee (apoapsis) and the position of a 
satellite. 
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APPENDIX A 


COMPUTATIONAL ALGORITHMS FOR GEOMETRIC LASER RANGING 


In this Appendix, the detailed computational algorithms are given for 
the solution of the multilateration equations and the sensitivity equations (See 
Paragraphs 7.2 and 7.3). For definiteness, the more general algorithms of the 
six-station configuration are given herein. Specialization to the 
four- and five-station cases or generalization to more than six-stations can be 
accomplished easily. The three-station collinear case is discussed in 
Appendix C. 


The input to the solution algorithm is the set of simultaneous slant 
ranges , i = 1 , 2, . . . 6; n = 1 , 2 , . . . N, where N is the total number of 
trajectory points used, and the nominal values of the station coordinates, namely, 

-y(o) y { ^ ) y(o) -y-(o) y(®) 7(0) -y (o) y(o) 7(0) -y(o) „(o) „(o)\ 

X 2 ’ X 3 ’ 3 ’ X 4 ’ Y 4 ’ 4 ’ X 5 ’ Y 5 ’ 5 ’ 6 ’ 6 ’ 6 )’ 

output will be the set of 12 relative station coordinates, (X^, X^, • 


. Z^) and, 


as a by-product, the satellite positions (x^, y , z n ), n = 1, 2, . . . N. The 
sensitivity analysis algorithm, on the other hand, requires as input the posi- 
tions of the satellite developed from an adopted orbit and the station coordinates, 
and calculates the magnification coefficients (sensitivity numbers) for the 
mapping of errors in the range measurements to the errors in each station 
coordinate. These algorithms are discussed below separately. 


A. 1 The Solution Algorithm 

In Paragraph 7.2, the basic idea of the method of solution, namely, the 
Newton-Raphson method applied to the multilateration problem has been 
outlined. Starting with the given ranges (i = 1, 2, . . . 6; n = 1, 2, . . . N), 


JPL Technical Memorandum 33-605 


A- 1 



and the nominal values of the station coordinates (X^, X^, 
the computation proceeds as follows: 

a) Calculate the range-square differences: 



). 


6. = p 2 _ p 2 ; , i = 2 , 3, . . . 6; n = 1, 2, . . . N. 

in in In 


b) Calculate the nominal satellite coordinates (n = 1, 2, . . . N): 



c) Calculate the following partial derivatives evaluated at the 

nominal values of the parameters. Note that the superscript 
for "nominal" is heretofore suppressed: 


8x X_ - x 

n c n 




dz 

n 

8X 2 
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d) Calculate the 3N x 12 matrix P(m, k) (m = 1, 2, 3, . . . 3N; 
k = 1, 2, . . . 12.) at the nominal values of the parameters, 
according to the following rules:* 


P(m, 1) 


9x 


- 2X . 


n 


j ax. 


ay 

- 2Y jax 


n 


- 2 Z 


9z 

n_ 

jax 2 


P(m, 2) = 


ay. 


dz 


7 v E. 7 7 E. 

" j ax 3 j ax 3 


P(m, 3) = 


ay. 


9z 


- 2Y jar;- 2Z j3Y. 


2 (X, - x ) if m < N 
4 n' 


P(m, 4) 


if m > N 


*The rules for setting up the [p] matrix follow from Eq. (7.2. 12). 


JPL Technical Memorandum 33-605 


A-3 



2(^4 - y n> 


i£ m - 


p(m> $) 


iim 


N 

if m > ^ 
< N 


4 «• 


p(m. 


if m 


> N 


7(Xc- x n 


if N < m 


i 2N 


p(m> 7 ) 


otherwise 


p(m> 8 ) 


\ if N - m 

[2 (* 5 " Yn* 

otherwise 


£ 2N 


_ „ x if N < m 

:z(Z r - z -' 


2N 


5 n 


p(mi 9) 


otherwise 

2N 


Um> 


“6 n 


p(m. 10 ^ 




p(m. 1^ 


otherwise 


if m > ^ 


otherwise 

> 2N 


in^u' z n ] ifm 


p(m. 


\- 


otherwise- 


jpL Technical 
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In the evaluation of the above matrix elements, the subscripts n 
and j appearing on the right-hand sides of the above equations 
are determined according to the value of m as follows: 


• 

If m < 

N, 

then j 

= 4, n 

= m 

• 

If N < 

m 

< 2N, 

then j = 

: 5, n = 

• 

If 2N < 

m 

< 3N, 

then j 

= 6, n 


m - N 


= m - 2N. 


e) Calculate the 3N x 1 column vector Q(m) at nominal values of 
the parameters according to the following formula: 


Q(m) = 2x X. + 2y Y. + 2z Z . - X 2 - Y 2 - Z 2 - 6- » 
n j 7 n j n j j j J jn 

where the subscripts j and n are determined as in d) above. 

f) Solve the linear system of equations expressed in matrix 
form: 


[P]AS = Q , 

for the 12 x 1 column vector AS by any standard technique. 

Note that for N > 4, the above system must be treated as an 
overdetermined system and the solution vector AS is to be 
interpreted as the least-squares solution. The components of 
AS are the corrections to the nominal station coordinates in 
first iteration. 

g) The corrected station coordinates S° + AS are treated as new 
nominal values and the algorithm returns to b) through f) for 
another iteration until two successive iterations agree to within 
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the specified accuracy. The final answer is the last-corrected 
solution. Note that the satellite coordinates (x , y , z ) are 
computed by the formulas of b) above, except, that the newly 
corrected values of X^, X^ and are used in place of their 
nominal values. 

The above algorithm was coded in Fortran. Computational 
experience indicates that the iteration converges rapidly. If the starting 
nominal values of the coordinates are accurate to within 10 meters, then only one 
iteration is required to generate a solution accurate to 1 cm. Even if the initial 
estimate is as poor as 1 km away from the true solution, only two iterations 
are needed to attain centimeter accuracy. 

A. 2 THE SENSITIVITY ANALYSIS ALGORITHM 

The basic ideas and the mathematical formalism of sensitivity, or 
error magnification analysis peculiar to the multilateration problem have 
been given in Paragraph 7. 3, where it was shown that the error dS in the station 
coordinate vector, S = (X^, X^, Y^, X^, Y^, Z^, . . . X^, Y^, Z^), due to 
errors in the measured range d is given by the solution of the linear system: 

A^dX,, + B^dX, + C^dY. + a ( ^dX. + b ( ^dY. + c^dZ. = dcr^ 
n2n3ninjnjnj n 

(j = 4, 5, 6; n = 1,2, . . . N), 

where the coefficients A^, etc., are defined by Eqs (7.3.4) through (7.3. 10). 

The computational algorithm for the simulation study outlined in Paragraph 7. 3 
proceeds as follows. 
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a) From the given configuration of stations and satellite 

trajectory points, calculate the slant ranges P- n and range 

square differences 6. by means of the formulas 

in ’ 


P 

in 


|(*„ - x 0 2 + (y n 




n 


1 , 2 , ... 6 ; 

1,2, . . . N, 


6. = P 2 - P Z . , i = 2, 3, . . . 6; n = 1, 2 . . . N. 

in in In 


b) Obtain a set of random errors dp. in the slant ranges from 

in 

a random number generator. These errors should be drawn 
from a Gaussian distribution with zero mean and unit standard 
deviation. The physical unit need not be specified for this 
purpose. 


c) Compute the coefficients 




•(j) a (j) 


n 


n 




4, 5, 6, 

1,2, . . . N, 


according to Eqs. (7. 3. 4) through (7. 3. 10), and the combined 
error vector daj^ according to (7. 3. 11). 

d) F orm the 3N x 12 matrix P according to Eq. (7.3.4) 
as follows: 
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e) Solve the linear system 


[P] dS = dQ 

by any standard technique to obtain the vector dS. The 12 

components of dS are the sensitivity numbers corresponding 

to the coordinates X_ , X,, Y_, X., Y., Z. . . . X,, Y,, Z, 

233444 obo 

respectively. 

f) Repeat the calculations b) through e) with a different set of 
random range errors a sufficient number of times such that 
meaningful statistical information, i.e. , the mean and the 
variance of each sensitivity number dS^> can be extracted. 

This algorithm has been coded in Fortran. Typical study results 
are discussed in Section V. 
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APPENDIX B 


DEGENERACIES IN MU LTILATERATION THEORY 

It has been repeatedly mentioned in the main body of this report that 
there exist station-satellite configurations such that the system of mathematical 
equations describing the multilateration technique is degenerate, i. e. , a unique 
solution to the station coordinates cannot be found. For example, when 4 or 5 
stations lie on a plane, then the system is always degenerate. This so-called 
coplanar degeneracy has been discus sed briefly in Paragraph 7. 4. In this Appen- 
dix, the nature of this degeneracy and why removal is possible using 6 stations 
is examined in more detail. In so doing, an exact solution to the coplanar 
6 -station problem will be derived. The conditions for further station degen- 
eracies due to peculiar station coordinate combinations also will be discussed. 

Another kind of degeneracy which arises as the result of certain 
geometric relationships between the trajectory points and/or the stations will 
then be considered. Degeneracies due to inadvertent use of badly correlated 
data are also briefly discussed. It should be pointed out that a complete, 
exhaustive treatment of all conceivable degeneracies has not been performed 
at the present time. However, it is believed that the most important degen- 
eracies have been isolated. Any further degeneracies that occur as the result 
of a peculiar interplay of configurational elements will show up as singularities 
in the sensitivity matrix or as nonconvergence of the Newton -Raph son iteration 
scheme (Paragraph 1 . 2 ). Finally, an operational technique for circumventing 
degeneracies will be explained. 
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B. 1 


Resume of the Coplanar Degeneracy 


In Paragraph 7.4, it has been shown that the fundamental equation of 
multilate ration, Eq. (7.2. 11), (also displayed in Figure 7.4. 1 for the first 
4 stations and 1 time-point), collapses for Zj(j=4, 5, ... I) equal to zero, to the 
following form; 


a . + (3. 
J J 



3n 



|n=l, 2, 
\3=4, 5, 



A(B - 1 ) 


where the quantities 6^ n are functions of the measurement data defined by 
Eq. (7.2.6), and the quantities a-, p-, and y- are functions of the station 

J J J 

coordinates given by the following expressions:' 1 ' 


“j * X j + Y j - X 2 X j - < X 3 + Y 3 - X 2 X 3> A(B - 2) 

p. = (X.Y 3 - X 3 Y.)/(X 2 Y 3 ) A.(B-3) 

y. = Y./Y_. A(B-4) 

J 1 3 

To recapitulate the argument presented in Paragraph 7.4 consider the following 
argument: suppose the system consists of only 4 stations all contained in a 

plane, then no matter how many measurements are taken, the system of equa- 
tions A.(B. 1) can at best determine three quantities**, namely the 3 station 


'Note that the values of X 2 and Y 3> which occur as divisors in these equations, 
cannot be zero, since they define the adopted coordinate system (See 
Section 4), 

Data point-caused degeneracies can further restrict the number of quantities 
determinable. See Paragraph B.4. 
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coordinate-parameter combinations a p^ and y^. These three quantities are 
not sufficient for the unique determination of the 5 unknown station parameters 
X X^, Y^, X^ and characterizing the planar 4-station configuration. Like- 
wise, in a 5-station coplanar system, at most 6 quantities are available for 
7 unknowns. Only when 6 stations are used, enough quantities (9) can generally 
be obtained for a unique determination of all of the 9 coordinate parameters. 
Singular configurations, however, could still exist even if there are 6 or more 
stations in the system. To discover such peculiarities, it is necessary to solve 
the coplanar 6- station problem exactly, as is accomplished in the next 
section. 

B. 2 EXACT SOLUTION OF THE COPLANAR SIX-STATION PROBLEM 
AND THE CONDITION FOR STATION DEGENERACY 

The previous section shows that in general nine quantities can be deter- 
mined from the planar multilate ration equations for 6 stations. These are the 
quantities a y p^, and y^, j=4, 5, 6 defined by Eqs. A(B-2), A.(B-3) and A.(B-4). 
These equations can be inverted easily. To proceed, consider the following 
operations . 

Equation A.(B-4) is rewritten as 

Y = y Y A.(B-5) 

J J ^ 

which upon substitution into Eq. A.(B -3) yields the following expression for X^ ; 

X. = p.X_ + v.X*. A(B-6) 

J 2 Y j 3 

Equations A.(B-5) and A.(B-6) are then substituted into A.(B-2), and after some 
algebra, the following result emerges; 
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p.(p. - 1) X 2 + Vj < Yj - 1) (X 2 + Y 2 ) + Zp. Yj X 2 X 3 = 


A(B-7) 


Equation A(B-7), written for j=4, 5, and 6, represents three linear equations in 

2 2 2 

the three coordinate combinations X + Y^, and X^X^. The condition for 
unique solution of these parameters to exist, which in turn determine X 2 ,X^ and 
is that the determinant 



P 4 (P 4 “ 

1) 

V Y 4 - 

1) 

2 P 4 y 4 

A s 

P 5 (P 5 - 

1) 

y 5 (y 5 - 

1) 

2 P 5 Y 5 


^6 - 

1) 

V V 6 - 

1) 

2P 6 y 6 


A(B-8) 


must not vanish. Note that and Yj are pure functions of the station coordi- 
nates defined by Eqs. A(B-3) and A(B-4). The same is true of the determinant 
The condition for station degeneracy ' 1 ' for the 6 station planar geometry can thus 
be stated simply as 


A(X 2 , x 3 y 3 , • 


•X 6 . y 6 > = °- 


A(B-9) 


Extensive numerical calculations have shown that degeneracies of this kind are 

infrequent. A general theory explaining the occurrence of station degeneracies 

has not been formulated at this time. 

To complete the solution of the problem when A 0, i. e. , when no 

2 2 2 

station degeneracy exists, the quantities ^X^, ^X^ + Y^, and £ = X -.X^ are 
determined by solving the linear simultaneous Eqs. A(2-7) written for j=4. 5, 
and 6 with the coefficients f3.(f3. - 1), etc., determined from the solution of 


"’'Degeneracy caused by the station configuration alone, independent from 
trajectory points. 
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A(B. 1) in terms of the measurement data. The station coordinates and 

are then given by: 


x 2 = ±JT 

X 3 = &//F A(B-10) 

Y 3 = ±(r\- t* 2 /i) l,Z . 

The ambiguity in sign in the above equations can be resolved once a definite 
coordinate system has been adopted. Having determined X^, X 3 and Y^ the 
coordinates of the fourth, fifth and the sixth stations then follow immediately 
from Eqs. A(B-5) and A(B _ 6). 

B. 3 ORBIT -DEPENDENT DEGENERACIES 

There exists the possibility that certain regularly related orbital 
trajectory points would cause degeneracy in the system equations. At the 
present time, no complete theory can be given for such degeneracies. However, 
certain obvious singularities can be isolated and should be avoided. 

The best way to discuss these orbit-dependent degeneracies is via the 
coefficient matrix P used in the sensitivity analysis described in Paragraph 7. 3. 
It was pointed out in Paragraph 7.3 that the system is degenerate if the coeffi - 
cient matrix is singular* . Inspection of this matrix, which has been explicitly 
written out in Appendix A, Paragraph A. 2, shows that it is singular if any one 
of the following submatrices P^ is singular, where 

'■'Singularity herein of a rectangular matrix is considered in the generalized 
sense of it being not of full rank . Thus, when a singular-value decomposition 
[5.3] i s performed on the rectangular matrix, less than the maximum allow- 
able number of singular values are different from zero. (The maximum 
allowable number of non zero singular values for a given rectangular matrix 
equals the smaller of the number of rows and the number of columns). 
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X 


1 


X. 

J 




P 


(j) = 


j=4, 5, 6. A(B-1 1) 


x. 


'N 


- X. 
J 




An obvious way for the matrices to be singular is for one of its columns to 
be a linear combination of the other two, or if constants kj and exist such 
that 


z - Z. = k.(x -X.) + k~(y - Y.) 
n j 1' n j 2 7 n j' 

for given j and all n. In other words P^ and hence P is singular, or the system 
is degenerate if all satellite trajectory points lie on a plane passing through one 
of the stations . This also includes the case in which all trajectory points lie on 
a straight line. This is true because a plane containing the straight line can 
always be constructed to pass through a given station. 

Other orbital degeneracies can conceivably occur. Computational experi- 
ence has indicated that using data from one pass of a circular orbit satellite 
generally leads to a singular P matrix. Even if multipass data (for the same 
circular satellite) are used, the sensitivity numbers obtained are rather poor. 
This indicates the near degeneracy of such configurations. However, if two 
circular orbits of different altitudes are used in conjunction with 6 stations, 
excellent sensitivities are generally obtained. The theoretical basis for such 
phenomena is not presently well-understood, and is proposed in this report as 
part of future study areas. (See Paragraph 10.2). 
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B. 4 


DATA DEGENERACIES 


Sometimes, due to the inadvertent use of highly correlated data, the 
resulting solutions to the multilate ration equations become ill-defined. For 
example, in obtaining the exact solution of the coplanar 6 -station problem in 
Paragraph A.2.Z, it is necessary to obtain uniquely the 9 quantities a\ , (3. , , 

j = 4, 5, 6. These quantities are determined by solving the system of equations 
A(B. 1). This process yields unique solutions only if the matrix D, where 


1 


1 


D = 1 




3 

'32 


A(B-12) 


6 2N 6 3N 


is nonsingular. Obviously, singularities of the D matrix are quite coincidental 
and the use of more data points (perhaps from another orbital pass) will 
generally remove the degeneracy. 

Since each data value is determined by the geometry of the station/ 
trajectory configuration, data degeneracies are merely different descriptions of 
orbit/station degeneracies discussed in the previous two sections of this 
Appendix. For instance, the matrix D above can be written with the aid of the 
multilate ration equations (7.2.3) and (7.2.4), as: 
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-2xiX 2 + X 2 
-2x 2 X 2 + X 2 


A(B-13) 


-2 Xi X 3 - 2 yi Y 3 + X 3 2 + Y 2 

-2x 2 X 3 - 2y 2 Y 3 + X 2 + Y 2 = D( v R.), 


where and R.. are the coordinate vectors of the trajectory points and the 
station locations. The condition that D be singular is thus equivalent to speci- 
fying certain station/orbit geometries. 


B. 5 OPERATIONAL METHOD FOR THE REMOVAL OF DEGENERACIES 
Despite the seemingly large number of degenerate possibilities which 
may cause system malfunctioning, there is, however, a practical method which 
will enable the analyst to discover these degeneracies and to avoid them. This 
simply consists of the combined use of the solution and sensitivity analysis 
algorithms detailed in Appendix A. To explain this, the reader is referred to 
Fig. B-l. Starting with Block 1, where the measurement data is assembled 
and a rough estimate of the station locating parameters in an adopted coordinate 
system is made, the Newton-Raphson solution algorithm is run (Block 2). If 
the algorithm converges to a solution, then the configuration corresponding to 
this solution is further checked by the sensitivity analysis algorithm (Block 3). 

If the analysis produces small error magnification factors, then an acceptable 
solution has been obtained (Block 4). Otherwise a nearly degenerate configura- 
tion is encountered (Block 5) and variation of the configurational element is 
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Fig* B- 1 . An Operational Plow Diagram for the Discovery 
and Avoidance of Degeneracies 
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indicated, ie. , the system must use data from other satellite positions and/or 
move some of the stations to other sites (Block 7). If however, the solution 
algorithm of Block 2 fails to give a convergent answer, then a degenerate situ- 
ation is obtained (Block 6) and the configuration should be varied (Block 7). 

The data and estimates associated with the new configuration are then of course 
reprocessed as indicated until an acceptable solution has been discovered. 
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APPENDIX C 


OBTAINING BASELINES USING ONLY 3 STATIONS 

In this Appendix, the following items relative to the three station 
solution are discussed in more detail. 

First, it is shown that if 3 stations can be aligned along a straight 
line, then simultaneous ranging to two satellite trajectory points will enable 
determination of the baselines to high accuracy. 

Second, analysis indicates that if the assumption of perfect alignment 
is relaxed, then an approximate solution to the 3 station (not perfectly aligned) 
problem, useful for surveying baselines up to 40 km, can be obtained. 

Third, simple formulas are given for the effect of small errors in 
alignment of the stations on the derived distances between them. It is shown 
that errors in the derived distances due to horizontal and vertical misalign- 
ments vanish if the airplane or satellite flies along certain optimal paths. 

Finally, it is concluded that the 3 station technique is best suited for 
measuring a network with mesh spacings of 10 to 40 kilometers, and for using 
an airplane as the retroreflector-bearing vehicle. 

C. 1 EXACT SOLUTION FOR PERFECTLY-ALIGNED STATIONS 

Suppose 3 stations are aligned in a straight line, as illustrated in 
Figure C-l and each of these stations ranges simultaneously to a satellite 
twice. Let a coordinate system with station 1 defining the origin, and the line 
joining the stations defining the X-axis, be adopted. The coordinates of sta- 
tions 2 and 3 are (X^, 0, 0) and (X^, 0, 0) respectively. Let the coordinates 

of the two satellite trajectory points be (x^, y^ , z^) and (x^, y^, z^). By defin- 

til til 

ing p to be the range from the i station to the n trajectory point, the 
following equations are evident: 
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Figure C-l. Three Station Collinear System 


SATELLITE 

ORBIT 


p 2 = (x - x .) 2 + y 2 + Z 2 , n = 1, 2; i = 1, 2, 3. A(C-1) 

'in n x 7 n n 

Note that = 0 in the adopted coordinate system. 

At first sight, it appears that Eq. A(C-1) defines only 6 equations for 8 
unknown quantities x^, y^, z^, x^, y^, z^, X^ and X^. However, the important 

quantities X^ and X^ can nevertheless be uniquely determined. By defining 
the usual range square differences (See Glossary): 


c 2 2 
6. = p. - p, ; 

m r in r ln 


i = 2, 3, 


n 


1 , 2 , 


A(C-2) 


it follows that 

6. = (x - X.) 2 - x 2 = - 2x X. + X 2 . A(C- 3) 

m n i n nil 

Equation A(G-3) determines four equations in four unknowns, and can be 
solved easily. The solutions for X^ and X^ are*. 
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* 6 22 ~ 6 21^ 6 21 5 32 " 6 22 6 31^ 

’32 " 6 31^ 6 22 " 6 21 " 6 32 + 6 31^ 


X 3 - 


(632 ~ 631)^21 6 32 " 6 22 6 31 ) 

22 ' 6 21 )(6 22 " 6 21 " 6 32 + 6 31 ) 


These solutions are exact under the assumption of perfect station 
alignment. A simple sensitivity analysis also shows that the solutions are 
stable with respect to slight changes in the measured ranges, with very low 
error magnifications. 

C. 2 APPROXIMATE SOLUTION FOR NONPERFECT ALIGNMENT 

If Station 3 does not lie on the straight line joining Stations 1 and 2 
(the X-axis), it can be assumed to be on the XY plane. Let be assigned to 
its Y coordinate. It will also be assumed that the values of Y^ and the Y- 
component, y^ of the satellite trajectory points are approximately known. The 
following ranging equations can be written: 

6 9 = - 2x X_ + X^ A(C- 6 ) 

2 n n 2 2 


6 , = - 2x X, - 2y Y- + xl + y\ 

3n n 3 7 n 3 3 3 


Equation A(C-7) can be rewritten as follows: 


6 ' = -2x X, + X" , 

3n n 3 3 


6 , = 6 -,+ 2y Y., - Y 

3n 3n 7 n 3 3 


These equations are exact, whether Y^ is small or large. 
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The important approximation that 6^ can be replaced by 6^ n is now 
made. This approximation is valid if is a small quantity; that is, if the 
stations are approximately aligned. Equations A( C- 6) and A(C-7)can now be 
solved exactly. The solutions are identical with Eqs. A(C-4) and A(C-5). with 
the quantities 6 , n = 1 > 2, in these equations replaced by the corresponding 

t 

quantities 6 as defined by Eq. A(C-9). 

C. 3 ERROR ANALYSIS 

From Equations A(C-4) and A(C-5) formulas for the errors in the 
derived quantities and X^ as functions of the alignment error Y^ can now be 
derived. Notice, from the symmetry of the geometry around the X-axis, that 
y and z are interchangeable in all equations; and therefore no generality is lost 
by representing the alignment error by Y^ only, since similar formulas for the 
effect of can be written by inspection. If station 3 is slightly displaced from 
the X-axis, then the appropriate values of 6 for insertion into the above formu- 
las will be & 2n and Eqs. A(C-6) and A(C-8); but the use of 6^ and 6^ n 

computed directly from the ranges will introduce an error. 

To the first order, 



d In X 2 


2 A ( 2 * * * 6 22 ~ 6 21 ) + 2 A * 5 21 6 32 " 6 22 6 31^ 2 a(6 32 ~ 6 31) 

^ 6 22 “ 6 21 ^ ^ 6 22 6 32 " 6 22 6 31 ^ ^ 6 32 " 6 31 ^ 


2 ^ b 22 " 6 21 ~ 5 32 + 5 31^ _ 

^22 ~ 6 21 " 6 32 + 6 31^ 


A(C-10) 


The differentials A&^ n are produced by the small displacement Y^ of 

the third station from exact alignment on the X-axis. The values of 6^ n> includ- 
ing the effect of the displacement, become: 
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- 6 21 = + 2 x 1 X 2- X 2 


- 6 22 + 2 x 2 X 2 ~ X 2 


■6 31 = +2x } X 3 + 2 Yl Y 3 - (X 2 + Y 2 ) 


-632 - +2x 2 X 3 + 2 ^2 Y 3 ‘ ^ X 3 + Y 3^‘ 


Hence, to the first order, 


2 A6 21 2“ u 22 


= 0, 


2 A6 31 = ‘ y l Y 3’ 


2 A6 32 -y 2 Y 3’ 


and also, 


(632 - 6 31 ) = 2(Xj - x 2 )X 3 + 2(y x - Y 2 ) Y 3 

[^ 6 22 6 21 * ‘ * 6 32 ‘ 6 31 ^] = 2 * X 1 " X 2^ X 2 " X 3 ^ " 2 ^ y l " y 2 

^ 6 21 6 32 " 6 22 6 31 ^ = 2 ^ X 1 " X 2 ^ X 2 X 3 ” X 3 X 2 ^ 

+ 4 ( Xl y 2 - x 2 Yl )X 2 Y 3 + 2 (y x - j 
whereupon, again to the first order, the expression for AX^/X^ 


AX 2 I 2(x l y 2 ' X 2 y l )X 2 + (y l ' y 2 )X 2 y l ‘ y 


X. 


2 (x x - x 2 )(X 2 X 3 - X 2 X 2 ) 2 (X 1 ' X 2 


y l - y 2 


2 (x x - x 2 )(X 2 - X 3 ) 


/ Y, 


> Y 3 

^ 2 )X 
becomes 

2 

;) X 3 
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Putting all terms in the above expression over the common denominator 


2(Xj - x 2 )X 3 (X 2 - X 3 ), 

the expression becomes 


AX 2 [(*jy 2 - ^yj) + x 3 ( Vl - y 2 )] Y 3 


X. 


<*1 - x 2 )X 3 (X 2 - X 3> 


A(C-ll) 


The formula for AX 3 /X 3 is virtually the same. Only the signs are changed 
in the first and fourth terms of the expression for d f n X, giving 


AX 3 [( Xl y 2 - x 2 yj) + X 2 ( Yl - y 2 )]Y 3 
X 3 ( Xj - x 2 )X 3 (X 2 - X 3 ) 


C-4 CONFIGURATIONS THAT MINIMIZE ALIGNMENT ERRORS 


A(C- 12) 


It is shown in this subsection that the trajectory points can be chosen in 
such a way that the errors in the computed distances X 2 and X 3> due to align- 
ment errors Y 3 and Z can be made to vanish. 

If the satellite or airplane flies a planar course defined by the equation 


y = A + Bx , 


then 


AX, = 


x 2 (a + bx 3 )y 3 

X 3 (X 2 " X 3> 


A(C- 1 3) 


AX, = 


(A + BX 2 )Y 3 
X 2“ X 3 


A(C- 14) 
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Three cases of interest may be considered: 

1. If A = B = 0, thenAX 2 = AX ? = 0; 

physically, this implies that if the vehicle flies over the X-axis, i. e. , 
in the XZ plane, then the effect of station misalignment' on both X 2 
and X^ vanishes. * 

2. IfA + BX 3 = 0, then AX 2 = 0 (but AX 3 t 0); 

this means that AX 2 = 0 if the vehicle flies in the plane defined by 

y = B(x - X 3 ), 

with B arbitrary. 

In other words, if the trajectory plane is parallel to the Z-axis and 
passes through Station No. 3, then the effect of alignment error on 
the coordinate X 2 vanishes. 

3. IfA + BX 2 = 0, then AX 3 = 0 (but AX 2 i 0); 

namely, if the trajectory plane is parallel to the Z-axis and 
passes through Station No. 2, then the effect of alignment error 
on the coordinate X 3 vanishes. 

Combining the above results, it is easy to prove the following: 

If the vehicle flies in a straight line course, approaching to or 
receding from Station 3, then the error induced upon the calcu- 
lated value of X 2 due to any small horizontal and vertical 
misalignment of the stations vanishes in the first order; 


The reason why the effects of the alignment errors vanish if the r etror eflector 
lies in the XZ plane can be seen intuitively from Figure C-2. The ratio of the 
range which would be measured if Station No. 3 were on the X-axis to the true 
range from the position displaced by Y 3> is cost) 1 . The variation of cos 4> from 
unity is of second order in Y^. 
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Figure C-2. Physical Interpretation of Alignment Error 

likewise, if the vehicle flies in a straight line to or from 
Station 2, then the error in vanishes. 

The above theorem has been tested numerically, and has been found to yield 
extremely accurate results (Typical values have been quoted in Sec VI). 

C. 5 CONCLUSIONS 

It has been shown in this Appendix that a range-only system for measuring 
inter station baselines is possible with only 3 stations, provided that these sta- 
tions are approximately aligned. The accuracy of the 3-station system can be 
greatly increased by optimally choosing the trajectories such that the effects of 
station alignment errors are minimized. 
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The requirement of approximate alignment of stations restricts the 
3-station system to short baseline applications. For baselines of less than 
40 km, an airplane can conveniently serve as the retroreflector -bearing vehi- 
cle. The maneuverability of the airplane can then be used to advantage for fly- 
ing the optimum trajectories. Further studies should be made of this promising 
technique. 
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APPENDIX D 


INVERSE TRIEATER AT ION: A METHOD OF STATION LOCATION 

The general technique of geometric multilateration developed in Section 7 
of this report requires no prior knowledge of any station coordinates, other than 
a very rough initial estimate for the purpose of starting the Newton -Raphson 
iteration. This self-starting property of the general technique is important 
whenever a new geographic area is to be surveyed. However, it is usually 
beneficial to establish new benchmarks (landmarks) in the general area where 
several stations already have been precisely located. For purposes of locating 
a new station, it is always possible to repeat the general multilateration calcula- 
tion simply by using the new station as one of a six-station or four- station 
configuration. However, with this procedure the coordinates of some of the 
established stations will be redundantly calculated, and thus waste much com- 
putational effort. In this Appendix, a simple algorithm, based on the well-known 
trilateration principle, is given whereby coordinates of new stations can be 
readily obtained once the coordinates of 3 master stations have been firmly 
established by using the general technique • This algorithm 
provides an auxiliary, efficient method for locating new benchmarks, i.e. , 
stations. 

D. 1 THE SOLUTION TO THE GENERAL TRILATERATION PROBLEM 

The trilateration principle, as will be explained below, is well-known in 
orbit determination. It has been successfully applied to the tracking of 
synchronous artificial satellites by the Goddard Space Flight Center 
However, as will be seen in the following section, it is equally valuable as a 
technique for locating new station coordinates. In the present section, the 
general theory of trilateration is discussed. 
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Expressed geometrically, the problem is as follows: to determine the 

three-dimensional coordinates of a point R(X, Y, Z) in space, given its distance 

p to 3 other points r = (x , y , z ) , n = 1 , 2 , 3 , whose coordinates are 
r n — n n n n 

assumed to be given in an arbitrary reference frame. 

The solution of the above problem is straightforward. Defining the 
vector slant - range _£ n by 

P = r - R , A(D-l) 

— n — n — 

the following equations may be obtained: 

p 2 = r 2 + R 2 - 2 r • R , n = 1, 2, 3. A(D-2) 

r n n — n — 

Equation A(D-2) represents a system of 3 equations in three unknowns, 
i.e. , the X, Y, Z components of R. The solution to these equations is exact 
and is given in algorithmic form as follows: 

(i) By use of the auxiliary quantities: 

x = x -x, y =y -y, z =z -z 
mn m n mn 1 m ' n mn m n 

c 2 2 , , 

6 =P - p. , m = 2, 3, 
m r m r l 

V = 2 l r m - r l - 6 rJ ’ m = 3 - 

D = X 21 y 31 ' X 31 y 21 
Pi = <“2 y 31 - "3 y 21 )/D 
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P 2 = (a 3 X 21 ' a 2 X 31 )/D 
q l = (z 21 y 3 1 ' Z 31 y 21 )/D 

q 2 = ^ z 3 1 X 21 " Z 21 X 31^ D 

, , 2 , 2 
a = 1 + qj + q 2 

b = " 2(pjqj + p 2 q 2 - Xjqj - y^ + Zj) 

2 2 2 2 
C = P 1 + p 2 ' 2 X 1 P 1 ‘ 2y l p 2 + r l ‘ P 1 * 

(ii) the solution is given by: 

Z = (- b ± J b 2 - 4ac) / (2a) 

x = Pj - q 2 z 
Y = p 2 - q, Z. 

The ambiguity in sign in the Z equation above must be resolved by appeal- 
ing to physical principles. As illustrated in Fig. D-l, it is evident geomet- 
rically that for a set of 3 given points in space and 3 ranges, there always 
exist, two trilateration solutions symmetrically located on both sides of the 
plane containing the 3 given points. Of these two solutions only one has 
physical meaning. For example, if the technique is used for orbit determination, 
then the physical solution is the solution that is above the plane of the stations 
(away from the Earth). 
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PLANE OF THE 3 
GIVEN POINTS 


Fig. D-l. Relation between the physical and the image solutions 
to the trilateration problem 


D. 2 INVERSE TRILATERATION FOR STATION LOCATION 

It is clear that the geometric statement of the trilateration problem given 
in the previous section makes no inherent distinction between the 4 points 
(3 known and 1 unknown) in space. By treating the 3 known points as master 
station locations, a trilateration will yield the coordinates of the fourth point 
(the satellite coordinates at a given time). Inversely, if the 3 known points 
are considered to be satellite trajectory points, then the fourth point, the 
location of an unknown station, can be trilaterated. 

The above trilateration scheme is especially simple when expressed 
in the adopted geometric coordinate system (Section 4.4), where the coordinates 
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of the 3 master stations are given as (0, 0, 0), (X^, 0, 0), and (X^, Y^, 0). 

Let the unknown coordinates in the master frame of reference of a new station 
be (X, Y, Z). Now since a sequence of satellite position coordinates (x n> y , z^) 
can be generated by direct trilateration from the master stations, then in the 
familiar notation of Section 7, these points can be expressed explicitly in terms 
of the known coordinates X^, X^ and Y^ and the range square differences 6^ 
(See Section 7. 2) as: 

= ( X 2- 0 /<2X 2> A(D - 3) 

V„ = ( X 3 + Y 3 - 6 3„ - ^„X 3 )/(3Y 3 ) A(D- 4) 

/2 2 2 \^^ 

z = + ( p. - x - y A(D-5) 

n \ In n } n / > ' ' 

where the positive sign of the square root of the z equation is chosen because 
the range point must lie above the XY-plane. Note that this solution, i. e. , 

Eqs. A(D-3) through A(D-5) is completely equivalent to the solution obtained 
by applying the trilateration algorithm given in the previous section (with the 
proper interpretation of symbols). 

Therefore, any three of the generated trajectory points can be used as 
the base for trilaterating the coordinates of a fourth point in space, i. e. , the 
new station location. As is evident, the entire algorithm of the previous section 
can be used for this computation. 

D. 3 POSSIBLE DIFFICULTIES 

The dual trilateration scheme explained in Paragraph D. 2 as in other geo- 
metric approaches, may run into difficulties because certain geometric configu- 
rations may result in singular algebraic systems. Therefore, care should be 
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exercised in the choice of the trajectory points used as the basis of an inverse 
trilateration. Generally speaking, a trilateration geometry is healthy unless 
all 4 space points lie on the same plane. Thus a satellite trajectory passing 
over a station should not be used to trilaterate the station location. Another 
singular situation is caused if the three trajectory points are collinear. All 
these difficulties can usually be overcome by data manipulation, A simple 
sensitivity analysis can assess the accuracy of any given configuration. This 
analysis has been performed, but is omitted herein. 
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APPENDIX E 


TIMING REQUIREMENTS 

The critical equation which relates the range to the total observed transit 
time is derived in this Appendix, i.e., Eq. (8.4.3) is explicitly obtained. In a 
stationary (nonrotating) geocentric coordinate system, let _r g (t) be defined as the 
instantaneous vector distance (magnitude r ) between the center of the Earth and 
the position of a satellite at time t. The instantaneous topocentric vector dis- 
tance from a ranging station to the satellite is denoted by p_(t) (magnitude p) . 
Finally, .R(t) is the instantaneous vector distance between the center of the Earth 
and the station on the surface of the Earth (magnitude R). 

Let a laser pulse be emitted at time t. It will arrive at the satellite at 
time t + A such that 


cA = |r g (t + A) - R(t)|. A(E-1) 

The pulse will now be reflected and arrive back at the station at time t + A + A' 
such that 

cA 1 = |r s (t + A) - R(t + A + A')|. A(E-2) 

Before solving Eqs. A(E-l) andA(E-2)it is important to know the magnitude 
of various quantities which are going to be used. The range will be of the order 
of 10 2 - 10^ km, so that the time A will be of the order of 3.10 2 - 30 • 10 2 sec. 
Hence, a Taylor series for r and R about t retaining only a very few terms will 
therefore be adequate. Accordingly: 

t' 2 

r (t + t') = r (t) + t 1 f (t) + r (t) + • • • 
s s s Z — s 

and 

' 2 

R(t + t') = R(t) + t' R(t) + R(t) + • • • 


A(E-3) 
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To obtain an estimate of the order of magnitude of the various terms in 
Eq. A(E-3) a realistic satellite orbit will be utilized. The radial velocity of a 
satellite is given by: 


whe r e : 



M 

G 


mass of Earth 
gravitational constant 


GM 


4 • 10^ [km^sec ^ ] 


A(E-4) 


r = r = radial distance from the center of the Earth, 

— s s 

a = semimajor axis of the orbital ellipse, 
e = orbital eccentricity. 

By choosing a representative orbit with minimum altitude above the 

surface of the Earth, H . = 900 km, and a maximum altitude of, H = 

min max 

45 00 km, it follows that: 

a = 9100 km e = 0.1978. A(E-5) 

From Eq. A(E-4)it follows that the maximum radial velocity (at r = 

s 

2 1 
all - e ]) is given by r =1. 325 km sec in which case r =0 and that 

7 max s 

the maximum radial acceleration (at perigee) is given by = -0. 726 • 

-3-2 

10 km sec in which case the velocity r g = 0. 

On the other hand the transit time of a pulse from station to satellite, 

_ 3 

when the satellite is at perigee, is about 15 • 10 sec, so that the satellite has 
_2 

moved 10 cm during this time. However, when the satellite experiences the 

largest radial velocity it can move during the then prevailing transit time of 
_3 

8 • 10 sec by 100 m; a change far from negligible. In any case range 
changes during pulse transit times are not affected by the satellite's accelera- 
tion within the accuracy desired (1 cm). 
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Stations located on the surface of the Earth will also move due to rotation 

-5 - 1 

of the Earth. Since the angular velocity of the Earth is w « 7. 3 • 10 sec 

an upper limit for both the absolute velocity and absolute acceleration is: 

• -1 •• -5 -2 

R = 0. 47 km sec and R = 3. 4 • 10 km sec respectively. As can be seen 
the acceleration again is totally negligible during the transit time of a pulse but 
due account has to be taken of the velocity. Neglecting the acceleration terms 
in Eq. A(E-3) it follows that 

cA = p(t) + r s (t) • p ( t) p (t) " 1 A A(E-6) 

cA' = P(t) + r g (t) • £(t)p(t) _1 A - R(t) • p (t) P (t) ~ 1 (A + A'), A(E-7) 

where of course 

p = r - R . A(E-8) 

— — s — 

It follows from Eqs. A(E-6), A(E-7), and A(E-8) that the observed total 
transit time A + A' = T is given by 

T(t) = | p(t)[l + c _1 p(t)]. A(E-9) 

Equation A(E- 9), the fundamental equation for the analysis of Section 8.4, shows 
that the total measured transit time T of a pulse emitted at time t is connected 
to the range at time t in the manner displayed. 
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APPENDIX F 


EXPECTED DATA ACCURACY 

As was discussed in Paragraphs 8.3 and 8.4 the process of obtaining 
synchronized range measurements will require that the range measurements or 
residuals be curve fit over suitable intervals by some assumed analytic func- 
tion. 

The question addressed in this appendix is related to the accuracy, 
i.e. , error bound, that will be obtained from the analytic model assumed for 
the curve fitting procedure when the number of measurements (range points) 
is large . 

F. 1 STANDARD DEVIATION OF DATA 

An important measure of the accuracy with which a given model, e. g. , 
polynomial will approximate the range history as a function of time, is the 
standard deviation of the raw data. If the modeled value of the slant range is 
denoted by p and the actual measurement is p , then the standard deviation of 
the model can be computed directly as: 

N 

V=2>n A(F - 1) 

n=l 

where N is the total number of range measurements taken in an interval of 
specified duration. 

F. 2 ASSUMED MODEL AND LEAST SQUARES FITTING PROCEDURE 

As was discussed in Section 8. 4, the data peculiar to this problem 
will be fit via least squares by the following function 

Ap = a + b(t - t ) + c(t - t ) 2 , n = 1, 2, 3, • • • , N A(F-Z) 

r n n o n o ' 
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where Ap^ are the true minus Keplerian range residuals, (See Glossary), 
th 

n refers to the n observation out of N total observations, t is the system 
time, and t is the adopted epoch time (See Para 8.4). 

As is well known from the theory of least squares, the object of 
obtaining the so-called best fit can be stated mathematically as that of minimiz- 
ing the function £2 defined by 


£2 


2 

* ( A n - a - b [‘i - l o] - - * 0 ] 2 ) 

+ (Ap 2 - a - b[t 2 - t D ] - c[t 2 - t D ] 2 ) + . . . A(F- 3) 

By extracting and equating the partial derivatives of £2 with respect to 
a, b and c to zero it follows that the least squares equations become 

N 

X^ Ap n 
1 

N 


N 

ill - hi 
1 

a - * of 

i 

N 

N 

N 

1 o 

i 

d 

L±i l 

W" 

XX ' l o] 2 
1 

E[‘n - ‘of 
1 

a - bf 

i 

ih - ‘oi 3 

i 

ih ■ *of 

i 


Xj A P n Dn " t o] 


1 

N 


X! ■ t o] 


A(F - 4) 


A useful simplification of these equations can be obtained by letting 

t = t + At(n - 1), 
n 1 

namely, by assuming that t is the laser data rate and choosing tj - t Q = 0, 
i. e. , the epoch and initial times are arbitrarily taken equal to zero so that 


t - t = At(n - 1). 
no ' 
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For large N (say, N > 50), the following approximations may be used: 




The above matrix equation can be solved and the quantities a, b, and 

c obtained as functions of N, Ap and At, i.e., 

n 

a = a(N, At, Ap), 

with similar equations in b and c. Furthermore the partial derivatives of the 
least squares coefficients with respect to the Ap-, which only appear linearly, 
can be obtained as 


3a 

9A Pn 


f(N, At), 


3b 

3Ap n 


= At) ’ w 


h(N, At). 


The algebraic manipulations become rather tedious and are omitted 
herein. Once the partial derivatives are obtained, the resulting error in the 
range assuming equal random errors in the Ap^ measurements, can be written 
as: 
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A(F-6) 


where cr^p^, the standard deviation of the data, is defined via Eq. A(F-l) with 
P n and "p n replaced by Ap^ and ^"P n * 

Note that the partial derivatives are summed from 1 to N with index n. 
Therefore using the previously stated summation formulas and after much 
patience 


3(3N - 7) 

4N(N 2 - 4) 


°Apf , 


A(F- 7) 


which for large N reduces to 


Ap 2 Apf 


A(F- 8) 


In summary, Eq. A(F-8) states that once the standard deviation of the 
actual data p n from the assumed curve fit "p n is computed, i.e. , once is 

determined, the actual value of the error in the range measurement can be 
obtained directly. It is this value, v&p’ that would be used in the explicit 
error mapping through the multilate ration equations to yield the random errors 


in the station coordinates. 
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APPENDIX G 


CURVE FITTING TECHNIQUE WITH RATIONAL FUNCTIONS 


As indicated in Section 8.3 there may be advantages in fitting residuals 
with rational functions. For purposes of outlining this technique only a trun- 
cated rational fit is discussed in this Appendix. 

If it is assumed that observations are related to time by means of the 
relationship 

p = A + Bt k , A(G-l) 

then the values of A, B, and k which make the square of any given residual, r, 
where 


r = A + Bt k - p , 


A(G-2) 


a minimum must be determined. In other words, for N observations, the 
function G, where 


£2 = r^ + 


+ r-, + 


+ r 


N 


A(G- 3 ) 


is to be minimized. 
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To accomplish this, the partial derivatives of f2 with respect to A, B 
and k must be set to zero. Proceeding formally, 






3A _ 

r l 

+ r 2 + .. 

• + r N = 

dn 


9r l ^ 

0r 2 ^ 

3B " 

r l 

3B 2 

SB + * * * 

dn 


9r l ^ 

3r 2 A 

3k 

r l 

8k + r 2 

3k + * * 


3r 


N 


'N 9B 
3r 


N 


'N 3k 


0 


0 . 


The expression for 3f2/3A implies that 

A + Bt^ - Pl + A + Bt^ - p 2 + ... = 0 



The expression for 3f2/3B implies that 

(a + Btj - Pj) ‘i + (a + Et^ - p 2 ) + . . . = o 

or 



A(G-4) 


A(G-5 ) 


G-2 
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Equations A(G-4) and A(G-5) form a linear system of two equations in A and B 
with the solutions 


A = 


( e *?) 2 - N ( E * r ) 


B = 


gqgj - 
x ^) 2 - N (^‘i k ) 


Hence, A and B are directly calculable once the value of k is known, 
expression for 3f^/3k implies that 

(A + Bt^ - P 1 )(Btj LOG^J + (a + Bt^ - P 2 )(Bt^ LOG e t) + ... 


or 


N 


N 


N 




LOG t. J A + 
e i 


^tfLOGeL B = ^ Pi tf 


LOG t. 
e 1 


Substituting for A and B from above yields 

( Ep fe ) . Ny &aLoG t 

(stf) -N(E‘f k ) 


(EtJiEpjtj) - (^i_ H £ p i l^ t k LOG t . ^ t k l OG t. 

' ei 11 ei 


A(G- 6) 
The 
= 0 

A(G- 7) 

A(G-8) 
= 0 . 
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The preceding equation is a non-linear equation in the parameter k. A 
solution can be obtained in one or two iterations via a numerical Newton pro- 
cedure using k = 1 as an initial estimate. 


G-4 
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APPENDIX H 


LASER PULSE SHAPE DEGRADATION BY MULTIPLE REFLECTION 
FROM CORNER REFLECTORS MOUNTED ON A SATELLITE 

Appendices H and I consider pulse shape degradation and pulse height 
degradation due to interference effects from multiple retroreflectors. It will 
be shown that for a large number of reflectors the interference effects will 
degrade the returned pulse, both in dynamic range (See Appendix I) and in 
pulse shape (Appendix H) to unacceptable levels. The analysis will indicate the 
need to use a small number of corner reflectors. 

In this Appendix the change of a rectangular laser pulse upon reflection 
from an array of retroreflectors will be analyzed via a simplified model. A 
numerical calculation for a typical set of parameters will also be presented. 

Suppose the corner reflectors are mounted uniformly and radially on a 
spherical satellite; i. e. , the reflectors are positioned equidistantly along a 
circle facing radially outward so that the angle subtended by the two lines 
emerging from the center of any two adjacent reflectors joining the center of 
the sperical satellite is a constant angle 0 q (See Figure H-l). 

14 - 1 

Next consider a ruby laser light pulse of frequency v = 4. 3 • 10 sec 
which is conveniently expressed analytically in an arbitrary coordinate system 
as : 


A(t) = [S(k • £ - wt) - S(k • £ - cot - T)] cos (k ■ £ - cot), A(H- 1 ) 
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where A is the amplitude of the electric field, 1c is the propagation vector, T is 

- 9 

the pulse duration (assumed equal to 2 • 10 sec), and S is the step function 
defined by: 

(1 if X > O 

S(X) = { A(H-Z) 

f 0 if X < O . 


The light pulse A(H-l), subject to the previously described reflector model, 
will therefore be reflected by the arrangement of mirrors. Assuming that the 
incoming light strikes one reflector normally, all the other reflectors will 
receive the light under a slant angle governed by 0 q . The intensity of the back 
scattered light is sensitive to the angle 0 between the axis of the reflector and 
the direction of the incoming light beam. 

It can be shown that, for reflectors currently in use, the amplitude of the 
reflected light may be expressed by: 


A = 



A(H-3) 


with a = 16°. The reflected light may now be expressed easily in a coordinate 
system centered at the center of the spherical satellite, viz: 


N 


A(t) = 2 + c COS nG o) - S ( t + t cos ne o - T )] 


n=0 


• cos (w— cos n0 - w t 
c o 


n 0 




\ ” a 

r 


A(H-4) 
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where r is the radius of the sphere and the summation index runs over all 
corner reflectors. 

Since the receiving photo multiplier is sensitive only to power, and since 
the inverse frequency v ^ is much smaller than the pulse length, a time aver- 
aging can be performed over the period of the carrier wave, i.e. , over the high 
frequency component. Hence, from Eq. A(H-4): 


2P(t) = 



( t + — cos n 0 
V c o 



r 

— cos n0 
c o 



e 



cos 



cos n 6 
o 


2 


+ {s [ S ( fc + C COS ne o) - S ( t + f cos ne o - T )] e 

. sin cos n ) | 


A(H- 5 ) 


where P(t) is the instantaneous power in the received signal. 

Figure H-2 displays a graph of Eq. A(H-5) for the following set of 

-9 15 -1 

parameters: T = 2 • 10 sec, r = lm, w = 2.84 • 10 sec and 0 q = 5°. 

Although only one set of parameters was used, albeit a typical one, it is clearly 
seen that the pulse shape degradation is significant. As is seen from Fig, H-2 
the leading edge of the reflected pulse is shifted somewhat in time. The con- 
stant fraction timing discriminator employed in the receiving electronics (See 
Section 9) integrates the pulse, determines the half width, and decides on the 
arrival time by selecting the particular point on the leading edge which corre- 
sponds to the half width. For the parameters chosen, it is clearly seen from 


H-4 
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Fig. H-2 that the arrival time determined in this manner is late by 
, -9 

1/8 X 10 sec; this corresponds to about 4 cm.' 1 ' This of course is unaccept- 
able for a system which purports to be accurate to the centimeter level. 

On the other hand, a shortening of the pulse would eventually create a 
separate return from each reflector, and therefore a series of individual but 
weaker pulses would be obtained. Again this is unacceptable since it would con- 
fuse the receiving electronics (See Section IX). Here, as well as in Appendix I, 
the necessity of a small number of corner reflectors is demonstrated. 


See A in Fig. H-2. 
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TIME, sec x 10 


Figure H-2. Pulse Shapes 
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APPENDIX I 


SIGNAL STRENGTH FLUCTUATIONS IN A LASER RANGING 
SYSTEM DUE TO OPTICAL INTERFERENCE BETWEEN THE 
MANY REFLECTORS ON A SATELLITE 


I. 1 INTRODUCTION 

An important parameter used in pulsed laser ranging systems is called 
the dynamic range . Leaving details aside for the present, the dynamic range 
in a sequence of M pulses emitted by a laser and returned to the receiving 
system^ is simply given by: 


dynamic range 


power of strongest returned pulse 
power of weakest returned pulse 


A(I-l) 


Because of interference effects, the pulses reflected from corner reflectors 
mounted on a satellite will vary in magnitude. In this Appendix the expected 
variation of pulse size will be calculated. The reason d'etre for performing 
this calculation is the simple fact that the receiving system, as it is conceived, 
consists of a photomultiplier (PMT) and a constant fraction timing discriminator 
(CFTD), and these components experience such variations . In this system the 
CFTD searches electronically for the complete pulse as generated by the PM. 
This search is performed in order to determine at which time the leading edge 
of the incoming pulse was received. In order to determine this time the dynamic 
range must be known because when the dynamic range is large, many of the 
received pulses cannot be handled by the receiving system, i. e. , they are lost 
for use in the process of range determination. 

In the following analysis the aim will be to determine the dynamic range 
for a reflecting system consisting of a large number of individual corner 

*The PM becomes saturated, i.e., overloads, at an input level of about 5V. 
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reflectors. The reason is twofold: first, the mathematical complexity is 

enormously simplified, and second, a large number of corner reflectors 
signifies a large return signal and therefore weaker laser pulses might be 
employed. The obvious advantages are: considerably reducing any safety 

hazards and also reducing the net cost per station. However it will turn out 
that the dynamic range is too large when the number of corner reflectors is 
large. Realizing that the dynamic range is just a measure of the interference 
effects caused by the interplay of the multitude of corner reflectors it is clear 
that a small number of corner reflectors will bring the dynamic range to 
acceptable levels. 

1.2 DERIVATION OF THE PROBABILITY DENSITY OF THE RECEIVED POWER 

The net signal received at the ground station will be the sum of the signals 
returned from the ensemble of retroreflector s. Since each retroreflector is 
at a slightly different distance from the ground station, the signals returned 
from the ensemble of retroreflectors will all have different phases. ' Conse- 
quently, the signals from the retroreflectors will constructively and destruc- 
tively interfere. 

In an actual system, the laser transmitter emits a sequence of pulses at 
a rate of 1 0 pulses /minute. Since each pulse is extremely short (<10 seconds), 
the relative positions of the retroreflectors (and consequently, the relative 
phases of the reflected signals) will not vary appreciably within the duration 
of a single pulse. However, the satellite will move a vast number of optical 
wavelengths between successive pulses. Consequently, the relative phases of 
the reflected signals will vary drastically from pulse to pulse, and the resulting 


14 


*At the light frequency (approximately 10 cycles /sec). 
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variation in interference effects will cause significant fluctuations in the net 
received power from pulse to pulse. 

1.3 ANALYSIS 

The linearly polarized electric field, E, of the light transmitted from the 
ground station can be represented as 


E T (t) = k cos (cot), 


A( I -2 ) 


where k is a gain constant and go is the light frequency. Now consider the light 

til 

reflected from the i reflector on the satellite. The E-vector of the light, 
when received at the ground station, can be represented as 


E.(t) 


a i E T 



A(t -3) 


th. 

where a^ is a gain constant, is the range between the station and the i 
reflector, and c is the speed of light. Combining Eqs. A(I-2) and A(I-3) 
yields 


E.(t) 


k a. cos co 
i 


(* 



A( I -4) 


*Only the magnitude (not the direction) of the E-vectors in these equations 
will be considered. 

**Throughout this analysis, the following effects are ignored: 1) atmospheric 

effects upon propagation, and 2) the variation of a^ with range. Since both of 
these effects will tend to increase the fluctuations of the received signal (as 
measured over a large number of pulses), the estimate of fluctuations will be 
conservative. 
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Since the argument of the cosine can be varied by any integral multiple 
of 2 it without changing the value of the cosine, the value of Zu>pjc can be 
replaced by 


4>, 



Jmod 2tt • 


A(I -5) 


Combining Eqs. A(I-4) and A(I -5) yields the final equation for the E-vector 

th. 

of the light which is reflected from the i reflector and received at the ground 
station: 


E.(t) = k a. coskot - 4>.) . 
1 l l 


A(I-6) 


since all E-vectors will add algebraically ’, if th'_ - a are N contributing reflec- 
tors on the satellite, the net E-vector of the light received at the ground station 
will be given by 


E R (t) 


E 

i - 1 


E.(t) 


•; - 1 


a. 

l 


cos (wt 


*i> 


A( I -7) 


*The aperture of the ground station's receiving telescope subtends a miniscule 
angle as seen from the satellite. Consequently, the light rays which are 
received from the ensemble of reflectors must be essentially parallel in order 
to be received. Therefore, all of the E-vectors will lie in the same plane 
i. e. , in the plane perpendicular to the direction of propagation, and the net 
received light will constitute a plane wave. If all reflected signals are linearly 
polarized the net received E- vector will be the algebraic sum of the individual 
E- vectors, and Eq. A(I-7) will be strictly valid. However, since mode-locked 
ruby lasers emit linearly-polarized light, and since all of the extant satellites 
employ metal-coated, e. g. , aluminized, reflectors which do not appreciably 
distort the polarization of the incident light, it is reasonable to assume that 
all reflected signals will be linearly polarized and that Eq. A(I-7) is valid. 


1-4 
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The assumption is now made that each reflected signal has the same 
amplitude (this will be essentially true if all reflectors have the same size and 
are mounted in the same plane), and therefore Eq. A(I-7) simplifies to 


E R (t) = ka 


E 

T - 1 


cos(ot - 4>j) . 


A(I -8) 


where each of the various a^'s have been replaced by the single parameter, a. 

In order to determine the probability densities of the amplitude and power 
of the net received signal, it will be necessary to specify the probability densities 
of the 4> i 's. First note that Eq. A(I-5) can be rewritten as: 


4 >- 


4it p i 


mod2ir 


A(I -9) 


where X. is the wavelength of the light. Since \ will have a typical value of about 
_5 o 

7 x 10 centimeters, i. e. , about 7000 A, even slight changes in (the range 
between the station and the i reflector) can cause large changes in 4irp^/X.. 
Since <J>. is equal to the value of 4-rrp./^- taken mod 2 tt, the various 4b' s will 

l r i l 

assume essentially- random values from pulse to pulse. In particular: 

a) the value of 4b during one pulse will be independent of the value of 
4>j during all preceding pulses (where j is any integer from 1 to N), 


and 


b) 


the values of all 4>^ 1 s will be independent of each other during any 
one pulse. 


*The relationship between p^ and 4^ is analogous to the relationship between x 
and y, where y is the 10th significant figure of log x. In general, even if a 
succession of x values has a definite pattern, the corresponding succession 
of y values will be essentially random. 
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Now consider the random process which consists of a large number of 
successive received pulses. In this process, each of the 1 s will represent 
a random variable, i. e. , each will assume random values from pulse to 
pulse. Since each <j>j is equally likely to have any value from 0 to 2ir, the 
probability density of each 4^ will be: 

P<*i> = SF = 0 s +i s 2 ” • 

Furthermore, since the various 4^ 1 s are independent, the joint probability 
density of the <f>. 1 s will equal the product of the individual probability densities, 
viz. , 


p(4>j * 4> 2 .... 4> N ) = pt+j) p(4> 2 ), • • . , p( 4> n ) • 


Equation A(l-8) shows that the net received E-vector is proportional to 
the sum of N sinusoids which have the same amplitude and frequency but differ- 
ent phases. The remaining problem consists in determining the probability 
distribution of the electric field Eq. A(I-8) and that of its amplitude. According 
to the foregoing assumptions the probability of E 0 (t) having the value E is given 

XV 

by the phase space available to all random variables <j>^ subject to the condition: ’ 

N 


E = 




cos(wt - <f>.) . 


A(I-10) 


i = l 


^Equation A(I-10) is conveniently normalized by the number N of corner 
reflectors to ease a subsequent passage to a large number of such reflectors. 
Note that E is a normalized electric field. 


1-6 
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Taking into account condition A(I-10), the probability density of E (t), can 

X\ 

then be written as 


P(E, t) = ( 2 tt) 


■ N 


2tt 

n 


d d> 6 
n 


N 


— / cos(u>t - 4>i) 

x /NX — i 


A(I-ll) 


where 6(x) is the Dirac 6-function. 

The integral A(I-ll) can be evaluated with the help of the integral 
representations of the 6-function and the Bessel function , yielding: 


P(E,t) 




A(I -12) 


*The 6 -functions in Eq. A(I-ll) may be thought of as the conditional probability 
density that Ej^(t) will have the value of E, given the values of <|>j , <t>2- • • ■ 
Hence, the probability that Ej^(t) will have any value is given by 



E - — Zcos (cot - 4> • ) 
N 1 


dE 


1, 


which is consistent with the definition of the 6-function. 


** These integral representations are as follows: 

/• oo 

6( y) = 27 / dxeixy 

J - oo 

and /» 2tt 

t / \ i / iy cose 

J o (y) 17 ' de e 


These formulas can be looked up, e. g. , in Ref. [ A. 1 . 3 .] 
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where is the Bessel function of zero order. Equation A(I-12) determines 
the instantaneous probability density of expression A(I-ll) and is seen to be 
time independent. The determination of the probability for the instantaneous 
electric field given in Eq. A(I-12) is not yet quite what is needed. The PM 

..r, 

tube responds to the average power' 1 ' W rather than to the instantaneous electric 
field. The derivation for the probability density of the average power W is 
somewhat more involved but follows strictly along the same lines as the deriva- 
tion given above. This is amply explained in [A. 1.2] and [a. I. l] . For the sake of 
brevity and for the simple reason that essentially no additional physical insight 
is gained from a straightforward but tedious derivation, the result will merely 
be quoted, namely, the probability density for the average power W is given by: 


J 0 ( Y n/2W) 


For a large number of corner reflectors, N, Eq. A(I-13), expanding the 
Bessel function, results in 

- 2 W 

P(W) - 2e . A(1 - 14) 


This equation shall be used in the subsequent analysis. 

1.4 PROBABILITY DISTRIBUTION OF THE DYNAMIC RANGE 

The probability distribution of the dynamic range of a sequence of M pulses 
will nowbe derived. 

Let the following definitions be introduced: 

v = power of the strongest pulse in a sequence of M pulses, 

^Averaged over the period of the carrier wave. 

❖❖Throughout this analysis, we will assume that M > 2. 
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u = power of the weakest pulse in a sequence of M pulses. 

gM<u, v) = joint probability density of u and v. 

Since the power of each pulse in the sequence has a probability density 
of P(W), then: 

P(u) du = probability that any one pulse will have a power between u 
and u + du, 

P(v) dv = probability that any one pulse will have a power between 
v and v + dv, 



probability that any one pulse will have a power between 
u and v. 


There are M possible choices for the weakest pulse, M-l possible choices 
for the strongest pulse (after the weakest pulse has been chosen), and M-2 
pulses whose powers have to lie between u and v in order for u and v to respec- 
tively correspond to the powers of the weakest and strongest pulses. Thus: 

M P(u) du = probability that the weakest pulse will have a power 
between u and u + du, 

(M-l)P(v) dv = probability that the strongest pulse will have a power 
between v and v + dv (after the weakest pulse has been 
chosen), 

-iM-2 

P(W) dW = probability that M-2 pulses have power between u and v 

J (assuming that the pulses are statistically independent). 
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Since the definition of g (u, v ) shows that, 

g (u, v) dudv = probability that, in a sequence of M pulses, the power of 
the weakest pulse lies between u and u + du, while the 
power of the strongest pulse lies between v and v + dv, 
and since the pulses are statistically independent, it therefore follows that: 


g^(u, dudv = t MP (u) du ] [ (M - 1 )P (v) dv] 


M-2 


P(W) dW 


L J u 


A(I- 15) 


Combining Eqs. A(l-14) and A{1-15) shows that, when the probability 
density of power in a single pulse is given by A(I-14), the joint probability 
density described by g^(u, v ) becomes: 


g M (u, v) = 4M(M - 1) exp j-2 [(M - 1) u + v] } { 1 - exp [-2(v - u)] } M 2 . 

A ( I -16) 


I. 5 Probability That the Dynamic Range of the Sequence of M Pulses will 
be 2 k. 

Let P^(k) be the probability that the dynamic range of a sequence of M 
pulses will be > k. Then, using the definition of g^(u, v) given above: 


P M (k) = Prob(v 2 ku) 


• oo /» oo 


g M (u, v ) dvdu > k - 1 


'0 J ku 


A(I -17) 


*Note that k must be 2 1 since we have assumed in the derivation of gj^(u, v) 
that v is 2 u. Specifically, it is assumed that v is the power of the strongest 
pulse and u is the power of the weakest pulse. 
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The integral A(I-17) can be evaluated and after some admittedly tedious 


algebra yields the result: 


P M (k) = 1 ' 


1 



M - 1 Terms 


A(I - 18) 


where A= (k- l)" 1 . Equation A(I-18) is displayed in Fig. 1-1 with M as 
a parameter. 

I. 6 CONCLUSIONS 

The receiving system can handle a ratio k = 50 efficiently. However, as 
seen from Fig. I - 1 , for a series of 100 pulses, the dynamic range k is twice 
as much with essentially unit probability. ''' This of course means that the pulse 
rejection rate is high. It must be remembered however that Fig. 1-1 is 
based on the simplifying premise of assuming a large number of corner reflec- 
tors. This assumption is tantamount to saying that interference effects are 
particularly severe. The use of a minimum number of corner reflectors is 
therefore a logical conclusion. 


til 

*This is just the dynamic range of Beethoven's 5 Symphony during the 
transition from the second to the third movement. 
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Fig. 1-1. Probability distribution of dynamic range 






APPENDIX J 


TRANSFORMATION OF COORDINATES 

Consider Figure J-I which illustrates the two basic coordinate 
systems of this study. In the Figure the numbers 1, 2 and 3 refer respectively 
to the three primary stations which define the geometric coordinate system. 



Fig. J-l. Inertial and Geometric Coordinate Systems 
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Figure 4, 5. 1 has displayed the basic relationships of the station 
geocentric coordinate vectors R^, namely [4. 7] . 


-I 


- - 

X 


cos 4> cos 0 

*■**/ 

Y 

= 

G^ cos <{> sin 0 

Z 

i 

G 2 sin 


where i is the station number and 


G ‘ ' [l -<2f - f 2 ) S in^J 1/2 +H 
(1 - f) 2 a e 

° 2 " [l -Uf-f^sin 2 *] 1 / 2 + H 


with 

= geodetic latitude, 

0 = X E + 0 g 

H = elevation of station measured normal to and relative to 

the adopted ellipsoid, 

X.£ = east longitude of station, 

0 s Greenwich sidereal time, 

§ 

a g = equatorial radius of Earth, 
f = geometric flattening of Earth. 

*v 

As can be seen from Figure J-l, is the vector from the dynamical origin 
to station one, i.e. , to the geometric origin. The vector from the dynamical 
origin to station 2 is denoted by R^. Subtracting these vectors and normaliz- 
ing results in 


/\ 


X 



(X 2 - X X )I + (Y 2 - Y^J + (Z 2 - Z l )K 
[(X^V+'V ? 1 ) 2 + '(Z 2 -Zj) 2 ] 1 / 2 


J-2 
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the overhead tilde denotes geocentric coordinates and where _I, _J, K are the 
unit vectors in the geocentric coordinate system, in this case, along the 
x-axis of the geometric system. Similarly, the vector subtraction of the 


third and origin station yields 


A 

s 


—3 ” —1 


~3 " — 1 


(X 3 - Xj)I + (Y 3 - Yj)I + (Z 3 - Zj)K 
[(X 3 - Xj ) 2 + (Y 3 - Yj ) 2 + (Z 3 - Zj ) 2 ] 1 / 2 


a unit vector toward station 3 in the geometric system. Taking the cross 
product of these vectors produces a unit vector in the geometric Z-direction: 

fi X3 = {[<y 2 - ?,)(z 3 - zp - (y 3 - yp(2 2 - zp]i 

+ [(Z 2 - Z 1 )(X 3 - Xj) - (Z 3 - Z 1 )(X 2 - Xj)]j 
+ [(X 2 - Xj)(Y 3 - Yj) - (X 3 - X 1 )(Y Z - Y 3 )] K | /| X X s| 


Hence using the definition of the cross product 

✓V /V /N . r /N 1 /? 

z = s x x s/[i -(x • sy\ ' 

where if all stations are in the northern hemisphere s = 1, if (X 2 - X^)(Y 3 - ) 


< X 3 - X ! )(Y 2 


Yj)_> 0; otherwise s = -1; however if all stations are in the 


southern hemisphere the previous rule is reversed unless station 1 is in the 


equator plane in which case _Z = R^/ 


R 


Finally taking the cross product of the two orthogonal unit vectors 


yields 


A 

Y 


✓V /\ 

Z X X, 


which completes the determination of the X, Y, Z basis. 

In terms of these unit vectors, the coordinates of the stations relative 
to the geometric system are now given simply by 
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X. 

= (R- 

-II)* 

X 

1 

— i 




r ^_ r 

/V 

Y. 

= (R. 


Y 

l 

— i 






Z. 

l 

= { -i 

" — 1 ) ' 

Z 


The transformation is now complete and can be expressed symbolically 


as: 


Polar Geodetic Rectangular Geodetic Adopted Rectangular Geometric 


4v 




H. 

1 




Y i> 


Z.. 

1 
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APPENDIX K 


MULTILATE RATION SYSTEM BIAS STUDY 


In this Appendix, the effect of possible station bias errors in range 
measurements on the resulting multilate ration solution is discussed. For this 
purpose, a simple model of the bias is used. It is shown numerically that 
arbitrary but small bias errors in the ranges measured at each station, that are 
maintained fixed throughout each individual calculation (but could be varied in 
different calculations), act very much like random range errors (of the Gaussian 
"noise" type) employed in the simulation results presented in Section V. A com- 
plete study of different kinds of biases will not be attempted here, but this study 
is recommended as a future task (See Section 10. 3, Task D). 

The reader is referred to Section 7.3 for the theory of accuracy simulation 
used in this report. Let the error in each measured slant range p^ n include a 
bias term b^ which is fixed for the i th stat ion, in addition to the Gaus-sian-noise 
type of random error gj n which has been used in the analysis of Sec 7. 3; namely; 


d P; 


in 


g. + b. 

' I VI I 


m 


A(K-l) 


i = 1 , 2, ... I, 
n = 1,2,... N. 


The "combined error vector" dcr^ VJ/ can now be separated into two terms as 
follows : , 


d <r 

n 


(j) 


y W + B W 

n r n 


A(K- 2 ) 
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j = 4, 5, ... 1 ; 
n = 1 , 2, . . . N , 


where 


(j) _ 


n 


F „ (r> s) P 


K fn g fn 


+ F j (£> §) 


P jn g jn 


I =1 


A(K- 3 ) 


P 


(j) 

n 


F „ §) 


p fn b e 


F. (r, S) p. b 
J - - K jn n 


= 1 


A(K-4) 


and the functions (r, S) are strict functions of the station- satellite geometry. 

The error dS^ induced on station coordinate S. due to ranging errors 
dp^ n is calculated via the matrix equation 


PdS = djr = Y_+3_, A(K- 5 ) 

where the matrix P is also a function of the geometry. Since Eq. A(K -5) is a 
linear matrix equation, we can separate the contributions to the induced station 
coordinate error into two parts, i. e. , "Gaussian noise contribution" dSy anc * 
"bias contribution" dS p : 


dS = dSy + dSp , 


A(K-6) 


where 


P dS y 

= Y 

A(K- 7) 

P dS p 

= P_ . 

A(K-8) 


K-2 
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It has been shown in Section V that random range errors (which make up 
the components of Y) resulted in excellent sensitivities in station locations for 
the 6- station configurations analyzed therein. It is reasonable to expect that 
bias errors of the type considered in this Appendix will behave similarly, 
inasmuch as the forms of the equations for the two kinds of errors are very 
much alike [Eqs A(K-3), A(K-4), A(K-7) and A(K-8) . This has been verified 
numerically. As a typical illustration, consider the 6-station configuration 
along the Southern California portion of the San Andreas Fault (See Para- 
graph 5. 1. 1). In addition to the 1 cm standard deviation Gaussian noise on the 
measured slant ranges, an arbitrary fixed bias (chosen also from a Gaussian 
distribution with zero mean and 1-cm standard deviation) was imposed on all 
the slant ranges measured from each given station. The induced errors on 
station coordinates are then calculated using the algorithm given in Paragraph A, 
Paragraph Al.2. The calculation was then repeated 10 times, each time with 
a different set of random noise and a new set of station biases. The results of 
these calculations are summarized in Table K-l. 

TABLE K-l 

The Average Error e and the Standard Deviation cr of the Coordinates of the 
San Andreas Fault Stations in Centimeters, Assuming 1 CM Standard 
Deviation Normal Random Errors in Ranging and 1 CM Standard 
Deviation Bias for Each Station 


Station 

€ X 

‘’’x 

£ y 

"y 

6 Z 

a z 

San Diego 

-0. 23 

2. 06 





Isabella 

-0. 50 

1.25 

0. 72 

3.94 



Santa Rosa I 

-0. 63 

1. 56 

-0. 21 

1.93 

-0. 83 

2. 03 

Millerton 

-0. 07 

1. 32 

0. 71 

4. 65 

0. 23 

2.21 

Blythe 

1 

o 

o 

2. 97 

1. 23 

5.03 

1. 64 

2.01 
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Comparing Table K-l with Table 5-2, it is evident that the accuracies 
of the computed station coordinates are not significantly degraded by the intro- 
duction of small station biases. In fact, the average errors remain well within 
1-2 cm, whereas the standard deviations (which measure the error for a single 
100-point calculation) are only slightly larger than the corresponding quantities 
in the absence of the biases. 
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APPENDIX L 


EARTHQUAKE HAZARD ESTIMATION ALONG THE 
SAN ANDREAS FAULT ZONE FROM A 
LONG- BASELINE STRAIN NET 

Clarence R. Allen, Don L. Anderson, and James H. Whitcomb 

Seismological Laboratory 
California Institute of Technology 

A major geophysical objective of accurately locating widely- spaced points 
as a function of time is to monitor the strain build-up in the earth's crust. Then, 
comparing the strain build-up with strain release resulting from earthquakes 
and from non- seismic slippage or creep, we can make a model of elastic 
deformation and its variation from place to place from which earthquake hazard 
can be estimated. Until now, geodetic techniques have not been able to measure 
relative locations over more than 100 kilometers to an accuracy of a few centi- 
meters, the accuracy desired for geophysical relevance. However, new tech- 
niques developed in the space program promise this precision. 

The San Andreas and associated fault zones in California and Mexico 
represent one of the most studied, least complicated (compaired to other fault 
zones) and potentially destructive seismic zones in the world. Therefore, this 
region provides a truly unique location for getting relative station motions and, 
thus, regional strain data which are both socially and scientifically pertinent. 

L. 1 Major Questions to Answer 

We must first determine the most important questions that can be answered 
by monitoring crustal locations in time and space. These questions should 
relate the California fault system to the new theories of plate tectonics and to 
possible earthquake prediction by strain monitoring. The questions are: 

(1) What are the ground movements at large distances (100-200 km) from 
the San Andreas fault zone? Knowledge of distant point movements is a basic 
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requirement in constructing a San Andreas zone strain model to be used as the 
basis of earthquake understanding. Current geodetic-net techniques, while very- 
accurate at short distances, are line- of- sight limited and are impractical beyond 
distances of about 100 km for the high accuracies needed. 

(2) How does strain build-up along the San Andreas zone vary as a func- 
tion of whether the fault is locked or unlocked; i. e. , how much of the accumulated 
strain is relieved by creep at the unlocked portions of the fault? It is now 
accepted that the type of slip that relieves strain in the San Andreas zone varies 
along the length of the fault. In some areas, much of the slip takes place as 
creep, that is, slow motion with little or no radiated seismic energy (earth- 
quakes). In other areas, almost all of the slip occurs rapidly in large but 
infrequent earthquakes. The former areas are considered unlocked inasmuch 
as the accumulating strain is being at least partially relieved as creep. The 
latter areas are considered locked inasmuch as most of the accumulated strain 
is being elastically stored for sudden release in further large earthquakes. 

Common sense dictates that the regional strain field along the San Andreas 
is strongly affected by whether the fault is locally locked or unlocked, and any 
strain measurement program should take this into account. For example, in 
areas of fault creep it is important to know if the input strain, measured by the 
displacement of points distant from the fault (100-200 km), is relieved by creep 
along or near the fault, which is relatively easy to measure. If the strain is 
relieved mainly by creep, then the likelihood of a large earthquake along that 
part of the fault is small. However, if only a small portion of the input strain is 
relieved by creep, then unlocked portions of the fault are also susceptible to 
large earthquakes. 

Figure L- 1 shows a map of the San Andreas fault zone and the regions 
of proven or suspected creep (unlocked) and regions with no observed creep 
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Figure L-l. Large Earthquake and Associated Faulting in the 
California and Nevada Region. The Locked and Unlocked 
Portions of the San Andreas Fault Zone are Delineated 
by the Shaded and Stippled Zones Respectively 
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(locked). Also shown are the larger earthquakes in the California region from 
Richter (1958, p. 472). Certainly if we take into account the fault breaks of the 
1857 and 1906 earthquakes (shown as solid lines in Figure A. 12. 1) the unlocked 
zone between San Francisco and Los Angeles has had fewer large earthquakes 
than the locked parts of the zone. This may be true of the unlocked zone of the 
Imperial Valley, although the picture there is not as clear. Whether this lack 
of seismicity is mainly due to creep relief of strain can be answered only by a 
direct monitoring of the input strain measured at distances far from the creep 
areas . 

(3) Is strain release occurring at places other than along the San Andreas 
fault zone? Two main possibilities exist. First, significant strain release may 
be occurring at other faults. As can be seen in Figure L-l, major earth- 
quakes have occurred on fault systems other than the San Andreas, and one 
should evaluate the importance of these events in terms of the regional strain 
build-up and release. Second, strain release may be taking place as permanent 
plastic deformation of the crust. This last statement, if true, contradicts a 
basic axiom of new plate tectonic theory, i. e. , that the plates are rigid. This 
axiom is by no means proven on a local scale such as in California and is an 
important point to investigate when developing a strain model. 

(4) What is the relationship of strain build-up along the San Andreas 
zone to sea-floor spreading in the Gulf of California? This question involves 
both the question of rigidity of the plates and the continuity of the Gulf of 
California's spreading rate; i.e. , is it spreading now at 6 cm/ yr as calculated 
from magnetic lineations in the oceanic crust? 

(5) What is the instantaneous velocity or displacement rate between 
adjacent plates? This is the most critical question that relates to earthquake 
frequency and hazard estimation along the plate boundaries. The displacement 
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rate determines how much slippage must be made up along the fault zone as 
creep or as earthquakes. The San Andreas fault zone forms the boundary 
between plates along its entire length, which explains its high seismicity. Other 
seismically active zones in the California-Nevada area also may be major plate 
boundaries, such as the fault system striking north through the Owens Valley, 
but data are incomplete for firm conclusions on these zones. 

L. 2 What is Needed? 

In order to measure a plate's horizontal movement completely, the time 

% 

history of the relative locations of at least two points in the plate is needed. 

Thus, with four points for two plates, a direct confirmation of the movement 
predicted by plate tectonic theory is possible. Only one point per plate is 
necessary if we assume some constraint on its motion, for example, that one 
boundary remain parallel to a strike- slip fault. 

The crustal motions that we anticipate are of the order of 6 cm/yr. 
Therefore, a measurement of the rate would begin to be significant if it could 
be made to within about 2 cm/yr. After the positions of the stations have been 
initially determined, the time to wait for a significant remeasurement depends 
on the accuracy of the measurement technique. 

Figure L-2 shows a schematic graph of a real station displacement 
rate r_ and two measurements of station position which will be within ±a and ±]d 
respectively of the actual position. The time interval between measurements 
is L Taking the worst cases from figure 1, we see that the rate r_ could be 
measured between 


a + b 

r = — • 

l 
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Figure L-2. Schematic Graph of a Real Station Displacement Rate r 
and Two Displacement Measurements with Errors ±a and ±b 
Respectively and With Intervening Time Interval i. 
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and 


r 


a + b 
i 


so that the error resulting from error a of the first measurement and error b 
of the second measurement would be 


a + b 
i 

Now assume that both a and b are 3 cm. The error in rate measurement is 
then 

± 6/i 

This means that we must wait three years before remeasurements would give 
an accuracy of ±2 cm/yr, but we will have increasing accuracy with each 
remeasurement thereafter. If we had to accept ±10 cm accuracy for both 
measurements then the waiting interval would be ten years to get the same 
rate accuracy of ±2 cm/yr. Current Jet Propulsion Laboratory work shows that 
±3 cm accuracy is now realistic. Therefore, a program begun now would 
measure the gross plate motions of the earth's crust to an accuracy of 2 cm/yr 
by at least 1976, with the accuracy improving at each successive measurement. 

L. 3 Recommended Station Placements 

In the framework of new plate tectonics, borders of active plates are 
mainly defined by zones of high seismicity, that is, high earthquake activity. 

This definition, when applied to the seismicity as represented by the largest 
earthquakes in the California-Nevada region in Figure L-l, determines 
three major plates in the area: the North American plate, the North Pacific 
Plate, and what we define as the Sierra plate as shown in Figure L-3. This figure 
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Figure L-3. Major Faults, Lithospheric Plates, and Proposed 
Stations for the Initial Long- Baseline Strain Net 
in Southern California. 
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also shows the major mapped faults in the area as dark lines. With this 
division we can readily see that a strain build up across the northern San Andreas 
zone, between the North Pacific and Sierra plates, can be different than that 
across the southern San Andreas zone, between the North Pacific and North 
American plates. 

Because Goldstone is to be one of the first stations due to its existing 
antenna facilities, it is best to have the first stations in the region which 
includes Goldstone, the southern San Andreas zone. Thus, our attention will 
be focused initially between only two of the above plates, the North American 
and North Pacific in southern California. We feel that the initial station loca- 
tions should be Goldstone, JPL, San Clemente Island, and the areas of Mount 
Laguna, and Glamis as shown in Figure L-3. The San Clemente- Mount Laguna- 
Glamis line crosses a San Andreas zone of known fault activity and suspected 
creep in the Imperial Valley of California. Indeed the Imperial Valley is the 
area of highest continuing seismic activity in the coterminous United States. 

Here, several parallel faults complicate the picture and the Mount Laguna - 
Glamis line straddles them for a measure of overall displacement across the 
valley. An added attraction in the Imperial Valley is that a unique history of 
geodetic data from government optical surveys is available as an independent 
check of the results. The San Clemente station to the west provides us with 
important data at a much greater distance from the fault zone. In addition, the 
relative displacements between San Clemente and Mount Laguna will shed light 
on current controversial questions of earthquake potential in the San Diego area. 

The Goldstone — JPL — San Clemente line is across a locked portion of the 
San Andreas zone, the Cajon Pass area, where many seismologists feel a great 
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earthquake (near magnitude 8.0) is possible within about 50 years. ^ In addition 
to monitoring strain build-up across this zone, it would be extremely important 
for us to measure these station locations before a large earthquake happens in 
order to make post- earthquake comparisons with theoretical strain models. 

When additional stations are possible, the next logical extension of the net 
would be north to the Sierra plate. Stations at Santa Ynez peak and north of 
Bakersfield shown in Figure L-3 measure displacements across the 1857 
earthquake break (solid line on Figure L- 1 ) of the San Andreas zone between 
the North American and Sierra plates. Motion of these two plates determines 
strain rates pertinent to the entire northern San Andreas zone including the 
San Francisco area. Further extension of the net with two additional pairs of 
stations across the northern San Andreas fault zone and two pairs south across 
the Gulf of California in Mexico would cover the entire San Andreas -Gulf of 
California plate boundary system. 

It is well to keep in mind in pioneer studies of this type that new and 
unexpected phenomena may make desirable a more detailed net in a critical area 
rather than an expansion to other regions of the San Andreas zone. Therefore, 
we should remain flexible as to future station locations. The portability of the 
proposed Jet Propulsion Laboratory system should make this not difficult. 


^For a review, see Brune, J. N. , 1971, Seismic Sources, Fault Plane Studies 
and Tectonics, Trans. Am. Geophys. Union, v. 52, no. 5, IUGG 178. 
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APPENDIX M 


PHYSICAL CONSTANTS 


Constant 

Adopted Value 

Earth Equatorial Radius, a g 

63781 65 meters 

Flattening, f 

1/298. 3 

Sidereal Rate, 0 

4. 37527 x 10 3 rad/min 

Speed of Light, c 

299714. 1 meters/sec 
, -6 

Second Harmonic, 

1082. 28 x 10 

Air Density (0°C, 760 mm), p 

a 

0. 001293 gm/cm 3 

Earth Gaussian Constant, k 

0. 07436574 (e. r. J 3 ^/ min 

Gravity Acceleration, g 

2 

980. 616 cm/sec 

Gas Constant, R 

8. 31 J/mole °K 
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APPENDIX N 


ABERRATION AND FRESNEL DIFFRACTION 

The need for a compact corner reflector has been demonstrated in 
Appendices H and I. Accordingly, the reflected light will consist of a narrow 
(~5 cm) diffraction limited beam. Furthermore, it is well known that the 
reflected beam will also suffer an angular shift due to the motion of the satel- 
lite. This angular shift is called aberration and its value, a, is: 

a = 2- , A(N-l) 

c 

where n is in radians and v is the velocity component of the satellite perpendic- 
ular to the line of sight. The factor 2 is due to the fact that reflected light 
rather than emitted light is being considered. 

Since satellites in circular orbits of from 500 to 1000 km have relative 
velocities approximately equal to 7. 5 km/sec, the maximum value of angular 
shift when ranging to these satellites is 

3 

| O' | = — = — * ^ ^ ^ — - = 50 pradians = 10 arcsec A(N-2) 

c 3 X 10 8 

This amount of angular shift is about 5 times the diffraction -limited beamwidth 
of a 5 cm corner reflector at ruby-laser wavelengths. Consequently if this size 
corner reflector is needed (in order to obtain sufficient capture area), the 
beamwidth of the corner reflector must be increased in order to cause a usable 
amount of reflected light to return to the ground station. 

One way of ameliorating the situation is to artificially widen the beam by 
means of a divergent lens or by rounding off appropriately the faces of the 
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reflector. The following analysis is a preliminary study of this problem 
using cylindrically symmetric apertures. 

The equation governing the diffraction pattern as depicted on Fig. (N-2) is 

derived in Sommerfeld's lectures on theoretical physics [*]. However, before 

going into a detailed analysis of the diffraction pattern corresponding to Fig. (N-2) 

a connection between Fig. (N-l) and (N-2) must be established. Figure (N-l) 

shows a corner reflector. A diverging lens has been placed in front of the 

corner reflector. The purpose of this lens is simply to make the reflected 

beam sufficiently divergent so that the receiving telescope, which incidentally 

is also the transmitting telescope, does not have to be aimed at a different 

angle upon receiving the returned pulse. To repeat, the aberration is an angu- 

_5 

lar deviation of about 5 • 10 radians assuming a satellite velocity of about 
7 km/sec. This angular deviation is sufficient so that a telescope transmitting 
a laser pulse will not be able to receive the pulse back unless the telescope is 
aimed at the new direction. This is very difficult to accomplish. Obviously, 
the alternative is to widen the reflected beam sufficiently to overcome the 
problem. 

In order to analyze this problem proceed from Fig. (N-l) to Fig. (N-2). 
Figure (N-2) shows (in cross section) a circular aperture illuminated by a point 
source at position S. The object here is to determine the far field diffraction 
pattern [R -*■ <» in Fig. (N-2)]. The light source at S (Fig. N-2) is just a replica 
of the virtual focal point of Fig. (N-l). The problem as stated in Fig. (N-l) can 
therefore be transposed to the configuration depicted in Fig. (N-2) and it is then 
amenable to mathematical analysis. 


Sommerfeld A. , Optik, Dieterich' sche Verlagsbuchhandlung, Germany, 
1950. 
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The analysis starts with Kirchoff' s Integral: 


U(P) 


■ * // 


dS 


ik(r+s ) 


r s 


[cos (n, r) - cos (n, s) 


A(N-3) 


Most symbols are amply explained by merely scrutinizing Fig. (N-2). The 
other quantities are defined as follows: A is the strength of the source [the focal 
point of Fig. (N-l)]. The wave length, X, of the radiation is assumed to be 
monochromatic. The angles (n, r) and (n, s) are clearly defined in Fig. (N-2). 

A number of assumptions are now made: 1) the radius of the circular aperture a 

is 2. 5 cm, 2) the distance r^ of the focal point from the aperture is large com- 
pared to its radius a, and 3) the distance R between the aperture and the obser- 
vation plane is essentially infinite. Therefore it follows that to a high degree of 
accuracy: 


cos (n, r) = 1 , cos (n, s) = -cos O'. A(N-4) 

Since the aperture of radius a is symmetrical about the z-axis it suffices 
to contemplate Fig. (N-2). Let 4 and q be the integration variables correspond- 
ing to x and y inside the aperture (z = 0). The location of S, the focus of the 
lens off the center axis, is 

S = (-x o , 0, -r Q ) A(N-5) 


The location of P is 


P = (tan o R, 0, R) 


A(N-6) 
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Therefore: 


= r 2 + ( £ + x ) 2 -f r) 2 
L o ^ o 1 


, 1/2 


= r 

J o 


r, 2 + 4 2 ^ X o 

+ ' + 

Zr r 

o o 


If x Q is not small, then: 


j.2 £,x 4x 6 

. o | o o 

r = r + + 


o 2r r _ 3 

o o 2r 

o 


x « r 
o o 


and 


s - JjR 2 + + (R tan a - l) 2 ^^ 2 - R (cos a) ^ - 4 sin a , 


where R — oo (far field). 
Let: 


c 2 , 2 .1 * , 
1+0 = 6 , 4 = 6 cos 9 


so that upon omitting all irrelevant phase factors: 


U(P) 


A 

2k 


1 + cos a 
cos 

Rr 

o 



ds exp 


| ~ik 4 sin 


a + ik 



+ 


x k 
o 

1 cos 

r 

o 


Letting 


lim A _ 

R — ♦ co R “ 


A(N-7) 


A(N-8) 


A(N-9) 



A(N-10) 
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and omitting all factors which vary slowly with o in the neighborhood of 

.4 

a « 10 radians, it follows that: 


U(P) 


r 

r 2 " 

<2 


/ x\ 

I 6 d 6 

1 d<j) exp 

ik -r 

c r 

exp 

ik 6 cos 4> (sin o- ) 

0 

Jo 

o 


\ r o/ 

i 


a(n-ii; 


Equation (N-ll) can be written immediately with the use of Bessel function 
representations as: 


U(P) 




A(N -12) 


Equation A(N-12) is true for large values of r Q (the focal length of the lens). 
Hence, it is seen that the diffraction pattern is just shifted by the amount x^/r^. 
This amount is a small number considering that is of the order of hundreds 
of meters. 

For the case in point x^/r^ therefore can be neglected. 

Figure (N-3) has been computed with values of the pertinent parameters, 
given by 

5 ------- - - - - - - 

ka = 2 • 10 , r = 250 m, 500 m, 1000 m, and ®. 

o 


The scale is normalized to the intensity P of the far field center line diffrac- 


tion pattern without a lens (r = co) and is plotted in db Idb = 10 log 


H) 


It is usual to call the case r = oo Fraunhofer diffraction since the 


o 

incoming light is a plane wave. The other cases deal with spherical waves and 
are called Fresnel diffraction; hence the caption on Fig. (N-3). Figure (N-3) 
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depicts the four cases just mentioned. Because of the rapid fluctuations and the 
wide spacing of the computed values, at least for the case r Q = it was not 
possible to actually draw a curve through the black dots representing the 
Fraunhofer case of Fig. (N-3). Indicated on the abscissa of Fig. (N-3) are the 
first three zeros of the r Q = oo pattern. It is seen that (at least for the aperture 
and wave length used here) the third zero appears in the neighborhood of 5 0 (jl 
radius. This is the critical angle for aberration, so that the beamcannot be 
seen by the receiving telescope. The other cases shown, although there are no 
nulls, are down in power by at least two orders of magnitude in the critical 
region around 50 pradians. This, however, constitutes sufficient intensity to 
be successfully processed by the receiver. 
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Figure N-2. The Mathematical Model Simulating Fig. N-l. Source is the 
virtual focus of the lens on Fig. N-l. The corner reflector is simulated 
by a circular aperture of radius a. The far field diffraction pattern 
is calculated by standard means (Ref. A(N-l). 
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Figure N-3. Frauenhofer (r = oo ) and Fresnel Diffraction. Explanations in the text 
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